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CHAPTER  1 

ELASTO-PLASTIC  CONSTITUTIVE  MODELS  AND  CONSISTENT  TANGENT 

OPERATORS 

1.1  Introduction 

The  study  concerns  the  nonlinear  representation  of  solids  through  finite  elements. 
In  such  an  analysis,  the  mathematical  description  of  the  material  behavior,  i.e.  the 
relationship  between  the  stress  and  strain  tensor,  a constitutive  relationship  is  required. 
Constitutive  models  (DiMaggio  and  Sandler,  1971,  Cap  Model,  which  will  follow) 
implemented  herein  are  of  the  elasto-plastic  variety  which  are  of  a phenomenological 
nature  and  formulated  from  experimental  observations.  It  is  not  realistic  to  try  to 
formulate  constitutive  models  which  fully  incorporate  all  the  interacting  mechanisms  of  a 
specific  material  because  any  constitutive  model  or  theory  is  a simplified  representation 
of  reality.  It  is  believed  that  more  insight  can  be  gained  by  tracing  the  entire  response  of  a 
structure  than  by  modeling  it  with  a highly  sophisticated  material  model  or  theory  which 
does  not  result  in  a converged  solution  close  to  the  failure  load.  An  important  objective 
in  the  present  study  is  thus  to  obtain  robust  numerical  tools,  capable  of  predicting  the 
behavior  of  the  structure  from  the  linear  elastic  range,  through  the  nonlinear  range 
including  the  complete  loading  history,  i.e.,  loss  of  strength.  To  achieve  this  objective  it 
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is  necessary  that  robust  constitutive  algorithms  be  complemented  with  advanced  solution 
procedures  that  result  from  the  finite  element  discretization. 

1 .2  History  of  Cap  Model  and  Solution  Procedures 

The  cap  model  originally  proposed  by  DiMaggio  and  Sandler  (1971)  consists  of  a 
failure  surface,  an  elliptic  hardening  surface,  and  a tension  cutoff  plotted  in  the  invariant 
stress  space.  The  cap  model  was  formulated  to  limit  the  dilatancy  which  developed  for  a 
number  of  elastic  perfectly  plastic  failure  surfaces  in  the  literature  (Drucker  and  Prager, 
1952).  The  proposed  cap  intersects  the  fixed  failure  surface  in  a non  smooth  fashion  and 
grows  or  shrinks,  i.e.  strain  hardens  or  softens  based  on  the  plastic  volumetric  strain  rate. 
Subsequently,  Sandler  and  Rubin  (1979)  modified  the  hardening  law  of  the  cap  to  prevent 
softening  response  which  would  occur  at  the  intersection  of  the  cap  and  failure  surfaces. 
Next  Simo  et  al.  (1988a),  using  the  concept  of  operator  split  (Ortiz  and  Pinsky,  1981; 

Simo  and  Taylor,  1986),  usually  referred  to  as  closest-point  projection  and 
unconditionally  stable,  formulated  the  scalar  equations  for  the  updated  stresses  for  all  of 
the  surfaces.  Expanding  on  this  work,  as  well  as  the  work  of  Simo  et  al.  (1988b), 
Hofstetter  et  al.  (1993)  developed  the  consistent  tangent  operator  for  all  the  surfaces  to 
preserve  the  quadratic  rate  of  convergence  in  a Newton  solution  procedure. 

Unfortunately  the  latter  development  requires  the  inversion  of  a 6x6  matrix  of  material 
stiffness  at  each  Gauss  point  per  iteration  which  could  prove  to  be  inefficient  and 
expensive  particularly  for  large  systems  (Boija  et  al.,  1991).  Hofstetter  recommends  that 
the  inversion  be  performed  algebraically  through  the  Sherman-Morrison  formula. 
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However,  since  all  of  the  components  of  the  material  stiffness  matrix  are  changing  at 
each  iteration,  the  latter  approach  has  no  direct  advantage  over  other  inversions,  i.e. 
multiplications  and  like  number  of  adds  (Press  et  al.,  1992). 

The  thrust  of  this  research  is  a direct  evaluation  of  the  consistent  tangential 
moduli  without  any  costly  material  stiffness  inversions.  The  moduli  are  developed  for 
each  surface:  failure,  elliptic  cap,  tension  cutoff,  as  well  as  comer  regions  by  determining 
the  active  yield  surfaces  from  Simo  et  al.  (1988b).  The  return  mapping  algorithms  for  the 
updated  stress  for  the  different  surfaces  are  found  from  the  closest  point  projection  (Simo 
et  al.,  1988a)  which  requires  the  solution  of  one  or  at  most  two  (elliptic  cap)  nonlinear 
scalar  equations.  Next  the  displacement  controlled  elasto-plastic  code  development, 
tangent  moduli,  and  quadratic  convergence  rate  are  presented. 

1 .3  Elasto-plastic  Finite  Element  Code 

Consider  a domain  acted  upon  by  external  forces,  such  as  body  forces  or  surface 
tractions.  For  the  body  to  be  in  equilibrium,  the  principle  of  virtual  work  states  the 
internal  work  must  equate  to  the  external  work.  The  basic  set  of  equations  to  be  solved  in 
terms  of  displacement  (d,+,)  are  (Boija  et  al.,  1989) 

^^IST  “ ^^EXr  )/it|  (1.1) 

where  (d„,.,)  is  the  internal  nodal  force  vector  and  (FExr)n+i  the  external  nodal  force 
vector  representing  the  integration  of  all  body  forces  and  surface  tractions  acting  on  the 
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body  for  the  n+1  load  step.  To  solve  Equation  1.1  for  the  unknown  displacements,  d,  the 
internal  nodal  force  vector  is  first  linearized  through  a Taylor  series  expansion: 


k*\. 
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where  d„.^|‘‘  is  the  nodal  displacement  vector  for  the  n+1  load  step  and  the  k"**  iteration. 
Knowing  the  internal  force  vector  to  be 

( Fjf^)  = fB„li  dV  (IJ) 


[B„+,  is  the  strain-displacement  matrix]  and  employing  the  chain  rule,  Equation  1.3 
becomes 


dF 
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dd. 


‘ fej.,  c,‘.  dv 


(1.4) 


where 


(1.5) 


is  referred  to  as  the  consistent  tangent  operator  obtained  by  evaluating  the  variation  of 
with  respect  to  the  strain  vector  Substituting  Equation  1.4  into  1.2  and 

setting  the  results  equal  to  Equation  1 . 1 gives  the  final  Newton  iterative  expression  for  d , 
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It  is  evident  from  Equation  1.6  that  the  system  of  equations  has  converged  and  the  body  is 
in  equilibrium  when  Fext  is  equal  to  Fnsn-  which  results  in  equaling  , i.e.,  the 
displacements  at  the  k+1  iteration  equal  the  displacements  at  the  k iteration  for  load  step 
n+1.  Usually,  however,  the  norm  of  the  incremental  displacements,  ||d'‘"^'  - d‘‘||,  or  the 
energy  norm  is  checked  against  a small  defined  tolerance  for  convergence. 

To  guarantee  the  quadratic  rate  (Equation  1 .2)  of  convergence  of  the  Newton 
solution  procedure,  the  consistent  tangent  operator  given  in  Equation  1 .5  has  to  be 
obtained  for  all  surfaces  and  comers  of  the  model. 

As  shown  in  Figure  1.1,  the  cap  model  consists  of  a failure  surface,  an  elliptic 
hardening  surface  which  is  a fimction  of  plastic  volumetric  strain,  and  a tension  cutoff 
surface  plotted  in  invariant  stress  space  (I,  and  / Jj).  The  first  invariant  of  the  stress 
tensor,  o,  is  I„  (I,=  o,(  + Oy  + ) the  sum  of  the  normal  stresses.  The  second  invariant, 

Jj,  of  the  deviatoric  stress  tensor,  s (o- 1, 1/3)  is  computed  from 


2 2 2 

xy  x2  yz 


(1.7) 


where  the  repeated  indices  (ij)  mean  summation  from  1 to  3 and  1 is  the  second  order 
unit  tensor.  Both  the  failure  and  tension  cutoff  surfaces  are  fixed  in  stress  space,  whereas 
the  elliptic  cap  expands  outward  with  loading.  McVay  et  al.  (1995)  reports  the  surfaces 
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Compressive 


Figure  1.1  Yield  Surface  of  the  Cap  Model 


and  the  associated  consistent  tangent  algorithm  for  Sandler,  Dimaggio  and  Baladi  cap 
model  as:- 


(a)  Failure  Surace  (non  hardening) 


/,(o)  = - a + yexp^^^'  - 0/, 


(1.8) 
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where  a is  the  vertical  intercept  of  the  failure  surface  on  the  second  invariant  Jj.  a defines 
the  shear  strength  of  the  material,  6 is  the  slope  of  the  failure  surface.  The  term  Yexp"***' 
is  used  to  provide  curvature  to  the  failure  surface.  The  consistent  tangent  for  the  failure 
surface  is 


da 
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(1.9) 


where^  - / Jj  n,.,  / ■/ J2.*Vi  ^nd  the  symbol  ® identifies  tensor  product, 
(b)  Tension  Cutoff  Surface  (non  hardening) 


/jCo)  = r-7, 


(1.10) 


where  T is  the  maximum  tensile  strength  of  the  material  obtained  from  experiments  and 
the  consistent  tangent  algorithm  is 
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(c)  Elliptic  Cap  Surface  (hardening  Surface) 


/2(o,e:0  = UrL? 


(1-12) 


X and  L are  the  hardening  parameters  and  locate  the  current  cap  surface.  X and  L are 
related  by  the  geometry  of  ellipse  so  that  only  X is  a true  hardening  parameter.  The 
horizontal  radius  is  (X-L)  and  the  vertical  ellipse  radius  is  (X-L)/R,  in  which  R is  a 
specified  parameter  defining  the  ratio  of  the  ellipse  radii,  if  R = 1 the  cap  surface  is  a 
circular  quadrant.  The  hardening  of  the  cap  is  controlled  through  accumulation  of 
compressive  volumetric  plastic  strains.  The  cap  is  located  by  specifying  X,  which 
establishes,  in  part,  the  initial  elastic  region.  In  turn,  the  elastic  region  provides  the  initial 
values  of  L.  The  consistent  tangent  algorithm  is 


(1.13) 


d)  Compression  and  Tension  Comer  Regions, 


(1.14) 
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In  Equations  1.7  - 1.14,  n is  the  unit  normal  tensor  to  s;  i.e.  n = s / V (s:s),  I is  the  4th 
order  identity  tensor  defined  as  1^^,  = ‘/2  (6jk  6j,  + 6j,  ),  6jj  is  the  Kroncker  delta  (2nd 

order  identity  tensor),  G is  the  elastic  shear  modulus  [G=E/2(l+v)],  and  K is  the  bulk 
modulus  [K=E/3(l-2v)].  The  coefficients  ai,Cj,anddj  involve  the  yield  surface 
geometry  (0,  y,  etc.)  and  the  converged  invariants  I,  and  J2„+,  (McVay  et  al.,  1995). 
The  detail  evaluation  of  the  consistent  tangent  operator.  Equation  1 .5,  along  with  the 
updated  stresses  based  on  the  closest  point  projection  are  presented  in  chapter  2. 

1.4  Organization 

Chapter  2 covers  the  formulation  and  implementation  of  an  efficient  time- 
dependent  viscoplastic  model  of  the  Perzyna  type  for  concrete,  rocks  and  soils.  The 
model  employed  is  an  inviscid  cap  model  first  proposed  by  DiMaggio  and  Sandler  (1971) 
and  later  by  Sandler  and  Rubin  (1979)  for  geomaterials.  The  model  did  an  excellent  job 
in  predicting  the  experimental  results  published  in  literature.  For  details  refer  to  chapter  2. 

Chapter  3 covers  the  behavior  of  adhesive  anchor  bolts.  Currently,  they  are 
designed  from  proprietary  tables  provided  by  the  manufacturers  based  on  the  laboratory 
pullout  tests.  Doerr  et  al.  (1989),  Cook  (1993)  and  Eligehausen  et  al.  (1984)  have 
developed  equations  to  predict  pullout  resistance  of  anchors.  These  approaches  require 
an  estimate  of  both  the  average  or  maximum  shear  stress  within  the  adhesive  bond  layer 
and  the  concrete  failure  cone  depths.  In  order  to  shed  more  light  on  the  behavior  of  these 
type  of  anchors,  analyses  were  performed  using  a state-of-the-art  elasto-plastic  finite 
element  program  and  compared  the  predicted  results  to  experimental  results.  The  study 
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revealed  that  uniform  bond  stress  applied  over  the  whole  anchor  did  an  excellent  job  of 
predicting  pullout  capacity.  Additional  studies  were  carried  out  for  different  concrete 
strengths,  anchor  diameters,  glue  thicknesses  , deeper  embedment  depths  and  different 
adhesive  strengths.  The  details  are  given  in  Chapter  3. 

Chapter  4 deals  with  a quadratic  axisymmetric  zero-thickness  element  of  the 
penalty  type.  Investigations  reveled  that  this  element  had  fundamental  kinematic 
deficiencies  inherent  in  its  response.  The  deficiencies  were  eliminated  and  an  improved 
quadratic  zero-thickness  was  formulated  and  implemented  into  the  nonlinear  finite 
element  program.  For  further  details  refer  to  Chapter  4. 

Chapter  5 covers  the  conclusions.  The  end  product  of  this  research  is  a PC  based, 
general  purpose,  nonlinear,  3-D  finite  element  program  which  can  handle  material 


nonlinearities. 


CHAPTER  2 

DEVELOPMENT  OF  A CONSISTENT  TANGENT  ALGORITHM  FOR  A 

MULTI  SURFACE  VISCOPLASTIC  MODEL  OF  THE  PERZYNA  TYPE 

2.1  Introduction 

The  thrust  of  this  chapter  is  formulation  and  implementation  of  an  efficient  time- 
dependent  viscoplastic  model  of  the  Perzyna  type  for  concrete,  rocks  and  soils.  The 
model  employed  is  an  outgrowth  of  the  two  invariant  inviscid  cap  model  first  proposed 
by  Dimaggio  and  Sandler  (1971)  and  later  by  Sandler  and  Rubin  (1979)  for  geomaterials. 
The  latter  has  been  noted  to  do  an  excellent  job  of  simulating  a number  of  features 
exhibited  by  many  geomaterials  (Simo  et  al.,  1988a);  however,  experimental  results  in 
concrete,  soils  and  rocks  indicate  that  such  materials  also  exhibit  significant  time 
dependent  behavior  even  for  small  loading  increments.  Failure  has  been  observed  to  be  a 
function  of  the  loading  rate,  i.e.  increasing  the  loading  rate  has  resulted  in  higher  failure 
stresses  and  diminished  failure  strains.  None  of  this  phenomeon  may  be  described 
through  an  inviscid  elasto-plastic  characterization. 

In  the  literature,  a variety  of  viscoplasticity  formulations  have  been  suggested 
(Katona,  1984;  Bodner  and  Parton,  1975;  and  Phillips  and  Wu,  1973);  out  of  these,  the 
two  most  popular  formulations  are  the  endochronic  theory  pioneered  by  Valanis  (1971) 
and  the  elastic-viscoplastic  theory  presented  by  Perzyna  (1966).  The  former  is  a 
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modification  of  classical  viscoelasticity  in  which  plastic  behavior  is  introduced  by 
intrinsic  time  scale  related  to  material  deformation.  Perzyna's  theory  is  an  outgrowth  of 
classical  plasticity  where  viscous  response  is  introduced  by  a time  rate  flow  rule  with  a 
plastic  yield  function.  Perzyna  viscoplasticity  theory  has  been  readily  adopted  because  of 
its  simulation  of  the  time  dependent  material  behavior  over  a wide  range  of  loading  and 
its  formation  is  readily  adaptable  to  numerical  algorithms,  suitable  for  the  finite  element 
procedure  (Katona,  1984). 

Some  of  the  soil  mechanics  specialists  who  have  worked  on  the  Perzyna  type 
viscoplasticity  are  Zienkiewicz  and  Cormeau  (1974),  Hughes  and  Taylor  (1978),  Pinsky 
et  al.  (1983),  Simo  and  Taylor  (1985),  Katona  and  Mulert  (1983),  Katona  (1984),  Desai 
and  Zhang  (1987),  and  Simo  et  al.  (1988b).  A brief  review  of  the  features  of  some  of 
these  works  is  presented  in  the  next  paragraph. 

Zienkiewicz  and  Cormeau  and  Desai  and  Zhang  implemented  the  Perzyna  type  of 
viscoplasticity  into  the  finite  element  program  based  on  an  explicit  integration  scheme. 

In  order  to  have  numerical  stability  very  small  time  steps  have  to  be  used,  which  makes 
the  algorithm  computationally  expensive.  Cormeau  (1975)  presented  time  step  limits  for 
stable  solutions  for  an  explicit  scheme.  Later  in  1978  Hughes  and  Taylor  offered  an 
implicit  scheme,  an  iterative  procedure  to  solve  nonlinear  equations  within  each  time  step 
to  eliminate  the  numerical  instability  due  to  the  explicit  scheme.  Katona  (1984)  proposed 
an  algorithm  for  the  Perzyna  type  viscoplasticity  model  which  involves  the  solution  of  a 
system  of  nonlinear  equations.  The  algorithm  relies  heavily  on  the  implicit  procedure 
developed  by  Hughes  and  Taylor  (1978).  Simo  et  al.  (1988b)  came  up  with  a new 
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algorithm,  based  on  the  closest  point  projection.  This  algorithm  results  in  a single 
nonlinear  scalar  equation  which  is  much  simpler  than  the  one  proposed  by  Katona 
(1984).  However,  the  Simo  et  al.  algorithm  has  never  been  implemented  into  the  general 
purpose  finite  element  programs  to  model  behavior  of  real  time  problems. 

This  paper  presents  an  algorithm  based  on  consistent  tangential  moduli  (without 
any  material  stiffness  inversions),  closest  point  projection  for  updated  stresses  which 
requires  one  or  at  most  two  nonlinear  scalar  equations.  The  algorithm  is  implemented  in 
general  purpose  finite  element  program  to  model  real  time  problems.  The  viscoplastic 
model  is  validated  by  taking  the  program  through  several  examples  mentioned  in  Katona 
(1984).  The  model  did  an  excellent  job  in  simulation,  was  unconditionally  stable,  and 
computationally  efficient. 

2.2  Formulation  of  Viscoplastic  Cap  Model 

The  elasto-viscoplastic  formulation  developed  within  the  context  of  the 
Perzyna  type  viscoplasticity  is  applied  to  the  Dimaggio  and  Sandler  (1971)  model.  It  is 
first  postulated  that  the  total  strain-rate  vector  € is  the  sum  of  elastic  and  viscoplastic 
components.  Thus,  we  have 

. . « . vp 

€ =€  +6  (2.1) 

where 

, e 

€ = elastic  strain-rate  tensor, 

. vp 

e = viscoplastic  strain-rate  tensor. 
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The  elastic  strain  rate  is  related  to  the  stress  rate  through 

o=Dee  (2.2) 

where 

D = elastic  constitutive  matrix 

Next,  the  viscoplastic  strain  rate  for  an  associative  flow  is  given  by  the  relation: 

(2J) 

where 

A,  = a material  parameter,  called  the  fluidity  parameter  (units  of  inverse 

time  and  denotes  the  relative  rate  of  viscoplastic  strain) 

f = viscoplastic  loading  fimction 
df 

„„  = current  direction  of  viscoplastic  strain  rate 

(j)(f)  = viscous  flow  function  (dimensionless  scalar  function) 

The  viscous  flow  function  <J)(f ) increases  monotonically  with  f and  defines  the 
current  magnitude  of  the  viscoplastic  strain  rate  and  should  satisfy  the  McAuley  bracket, 
which  is  given  by 

((J)(f))=  (j)(f)  for  f>0 

((J)(f))  = 0 for  f^  0 (2.4) 

The  next  most  important  thing  is  the  choice  of  the  viscous  flow  function  4>(f).  For  the 
two-invariant  constitutive  model,  two  most  commonly  used  forms  will  suffice  for  many 
geological  materials  (Katona  1984); 


WH^f 

Jq 
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(2.5) 


and 


(j)(/)=exp(^)^  -1 
Jo 


(2.6) 


where 

N=  material  parameter  (real  and  positive  number) 

fo  = normalizing  constant  (same  units  as  loading  function,  f ) 

2J  ViscQDlasticitv  Loading  Function 

The  viscoplasticity  loading  function  or  the  yield  function  for  the  Sandler  et  al. 
(1976)  model  consists  of  a failure  envelope  f,  (o),  a tension  cutoff  region  fj  (o),  and  a 
strain  hardening  cap  fj  (o,  e/" ),  Figure  1.1.  Both  the  tension  cutoff  and  the  failure 
envelope  are  stationary  surfaces  (ideal  plasticity),  whereas  the  cap  translates  based  on  the 
magnitude  of  the  viscoplastic  volumetric  strain  e/P  generated.  The  functional  forms  for 
the  surfaces  are 


/,  (o)  = ) for  T < I^<  L{ef)  (2.7) 


where 
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F,Ux)  = a - Yexp  *^^'  + 0/, 


(2.8) 


for  I(€7)^/,^X(e:^ 


(2.9) 


fia)  = r - /,  = 0 , for  f = T (2.10) 


where  a,  p,  y,  0,  and  R are  material  parameters  for  the  cap  and  failure  envelope;  T is  a 
material  constant  (maximum  tensile  stress  obtained  from  split  tension  test)  referred  to  as 
the  tension  cutoff;  I,  is  the  first  invariant  of  the  stress  tensor  o , tr  (o);  s is  the  deviatoric 
stress  tensor  defined  as 


s = a 


(2.11) 


1 is  a second  order  unit  tensor;  and  f Jj  is  the  second  invariant  of  the  deviatoric  stress 
tensor,  s.  As  generally  assumed  in  soil  mechanics,  positive  normal  stresses  are 
compressive.  The  point  of  intersection  of  the  cap  and  failure  surface  is  given  by 
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Z(e:0  - X(eT)  ^ RF^iL) 


(2.12) 


The  growth  of  the  elliptic  cap  is  controlled  through  the  viscoplastic  volumetric  strain, 
as  suggested  by  Sandler  and  Rubin  (1979), 


in  which  W and  D are  cap  hardening  parameters  which  are  important  for  capturing  the 
shape  of  the  stress  / strain  loading  curve.  W represents  the  maximum  volumetric  strain 
that  can  be  produced.  With  regard  to  D,  it  is  generally  recommended  to  choose  D Ix  I 
should  be  less  than  0.5.  Decreasing  D or  W results  in  increased  stress  magnitude,  i.e.,  a 
steeper  stress  / strain  loading  curve.  For  any  problem  these  values  can  be  obtained  by 
trial  and  error.  is  the  first  invariant  of  the  viscoplastic  strain,  e''’’ , tensor,  or  tr(  e'’’’ ). 
The  model  employed  herein  assumes  that  the  cap  only  hardens  for  stress  loading  onto  the 
cap;  no  softening  or  shrinkage  of  the  cap  due  to  loading  onto  the  other  surfaces  or  comers 
is  allowed.  Sandler  and  Rubin  (1979)  proposed  such  a characterization  for  geological 
materials,  and  rock  in  particular. 

During  viscoplastic  loading,  the  current  stress  point  (the  invariants  I,  and  / Jj 
which  are  used  to  define  the  current  static  yield  surface)  is  allowed  to  be  outside  the  static 
yield  surface  f = 0.  The  static  yield  surface  forms  the  boundary  between  elastic  domain 
for  which  f < 0 and  the  viscoplastic  domain  for  which  f > 0.  Viscoplastic  flow  will  occur 


W [\  - exp  ] 


(2.13) 


and  continue  to  occur  at  a constant  rate  if  the  yield  function  is  a nonhardening  nature,  and 
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on  the  other  hand  the  viscoplastic  flow  will  occur  at  a decreasing  rate  (such  that  fj(o,  e/’’) 
tends  to  zero,  which  means  tends  to  zero)  if  the  yield  function  is  hardening  in 
nature.  This  means,  the  static  yield  surfaces  move  out  on  a real  time  basis  to  eventually 
form  a static  yield  surface  containing  the  imposed  stress  state;  at  this  stage  we  say  the 
solution  is  steady  state  = 0.  The  resulting  strains  accumulated  during  this  loading 
would  be  identical  to  the  corresponding  inviscid  plasticity  solution. 


2.4  Consistent  Tangent  Operator  and  Updated  Stresses 

From  the  converged  solution  at  time  t = t„,  the  elastic  trial  stresses,  are 
computed 


= o„  . ; 


(2.14) 


where  is  the  elasticity  tensor  and  o„  is  the  converged  stresses  at  time  t = t„.  Ae^,  is 
the  total  incremental  strain  tensor  at  time  t = t„+,  . Knowing 


«...  = «.  " . C^:(Ae„.,  - AeJ.,) 


(2.15) 


and  making  use  of  equation  2.14,  equation  2.15  reduces  to 


«...  = C.  - C^:AeZ, 


(2.16) 


where  C^:A6„„|'^  is  the  plastic  corrector.  Substituting  in  for  , o^.  becomes 


«...  = «:..  - KAeZl  - 2GAe^ 


(2.17) 
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K and  G are  the  elastic  bulk  and  shear  moduli,  respectively.  From  equation  2.3, 
we  get  incremental  and  deviatoric  viscoplastic  strains  as 


Ae7  = 3 AA  ((})(/)> 


(2.18a) 


Ae 


vp 
/l  + l 


- AA<4>(/)) 


(2.18b) 


From  Equation  2.17,  stresses  at  time  t^,  reduce  to 


'n*\ 


= o 


n*\ 


AA((j)(/)) 


dfi 

3K—\ 

dl, 


+ 2G 


n 


/2 


(2.19) 


The  next  step  is  to  find  the  active  portion  of  the  yield  surface  for  the  particular 
load/time  increment.  The  active  portion  of  the  yield  surface  may  be  characterized  solely 
on  the  basis  of  the  elastic  trial  stresses  (Simo  et  al.,  1988b). 


.tr 


ntr 


I "■ 
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(2.20) 


where  / and  I|,n+i"  are  the  trial  invariants  which  are  computed  from  the  trial  elastic 
stresses.  The  value  of  1,  n+i*"  is  determined  from  the  tr(o„+|'" ),  and  -f  J2n+i“  is  found  from 


(2.21) 


where  the  deviatoric  stress,  s'" , is  defined  in  Equation  2.11.  Next  the  invariants  I, 
and  / J2  n+i  are  computed  as 


c,  - 9KAXm) 


dl, 


Jn*l 


- GAkim 


dfx 


■ n + l 


(2.22) 


Next  the  derivatives  of  fj  with  respect  to  the  invariants  are  evaluated  based  on  the 
active  portion  of  the  yield  surface.  The  number  of  unknowns  are  reduced  by  substituting 
one  in  terms  of  the  other.  Development  of  a consistent  tangent  operator  for  each  of  the 


five  plastic  loading  modes  (failure,  elliptic  cap,  tension  cutoff,  tensile  comer  region  and 
compressive  comer  region)  is  shown  below. 


Failure  envelope:  The  yield  surface  for  the  failure  envelope  is  characterized  by 
f\  = ^^2  ~ ^ ^ yexp  ^^^'  - 07,  (2.23) 

where  I , and  VJ  2 are  given  by  Equation  2.22.  By  taking  the  derivatives  of  fj  with 
respect  to  invariants  and  substituting  back  in  Equation  2.22, 

'..I  - Cl  ♦ 9/^AX<(K4)>(YPe  '’''-+e)  (2.24a) 

A...I  = '^•C  - GAA.<(ti(/,))  (2.24b) 


where  ((j)(f,))  is  given  by  equation  2.3 


(4>(/i)) 


/J2  - a + yexp'^^'  - 67, 
fo 


N 


(2.25) 


From  Equation  2.24a 


AA(4)(/;)) 
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Substituting  this  into  (2.24b)  results  in 


/J. 


2.W  + 1 


= /J'" 


9K 


id  -r  )] 

^ (Y  p e +0), 


(2.26) 


Once  / J2  „+i  is  expressed  in  terms  of  I|  , the  nonlinear  scalar  equation  2.24a 
for  I , „+|  may  be  solved  by  a Newton  iteration  technique.  / „+,  and  (4)(f|))  are  solved 
in  the  process  of  solving  for  I,  . 

The  consistent  tangent  operator,  Equation  1.5,  is  subsequently  obtained  by  first 
expressing  the  stress  tensor,  o^,  in  terms  of  its  invariants. 


o 


n*\ 


1 * n 


(2.27) 


and  then  differentiating  it  with  respect  to  e„+, , using  the  chain  rule 


da 


de 


= ll® 


dl. 


n*\ 


n*\ 


de 


V2/J. 


2,n*l 


n*\ 


dn 
'de~ 


r d/j,. 

+ /2  n<s 


n*\ 


de 


(2.28) 


n*l 


where  the  symbol  ® identifies  a tensor  product  in  the  sense  that  (a  ® b)ijy  = ajj  b^, , for  any 
second  order  tensors  a , b . The  partial  derivatives  of  the  invariants  with  respect  to  e„+, 
may  be  determined  from  Equation  2.24, 


3*:i  . 9/(AA.[e.vP«''‘''"“. 


de 


n*\ 


de 


rt+i 


1 + 9/:  AA.  (({)(/;))  YP^e'’^^'"*') 
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and 


= /2G„  - 


de 


n*l 


de 


(2.29) 


n*l 


In  Equation  2.29  the  derivative  3((l)(f|))  /de„^,  can  be  determined  from  Equation  2.25; 


-r; " 7 

^^n*l  f\ 


diJ. 


de 


de 


>1+1 


(2.30) 


By  substituting  Equation  2.29  into  Equation  2.30  and  rearranging  the  terms  we  get  an 
expression  for  3((J)(f,))/0e^,.  Plugging  this  expression  into  Equation  2.29,  the  final  forms 
of  ai,  /ae^,  and  a/Jj  lde„^^  are  given  by 


a/y. 


2.W+1 


de 


Cf,\^(C^^C,C^)n 


n*\ 


(231) 


where 


3*K 

^2 


Cj=v/2  B G — — 
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r- 

' {D,^\) 


D,(^+0) 


and 


A=-i^e 

B=9Kmm) 


^2=1 


D^=G^X{m)J 

J\ 


To  complete  the  evaluation  of  the  consistent  tangent  operator,  in  Equation  2.28, 


the  derivative  dn  /0e„+,  given  by  Simo  and  Taylor  (1985)  is  used: 
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dn  _ 


[/-  — I ® 1 -n<s>n] 


1?  3 


(232) 


where  I is  the  4th  order  identity  tensor  defined  as  = 1/2  (6;^  6^,  + ) and  6jj  is  the 

kroncker  delta,  6,^  = . Substituting  Eqs.  2.3 1 and  2.32  into  Equation  2.28  results  in  the 

final  form  of  the  consistent  tangent  operator  for  loading  onto  the  failure  surface, 


where^=A.n.i/A,\.i  • 

Elliptic  cap:  The  cap  is  the  only  portion  of  the  yield  surface  which  hardens  or 
moves.  In  order  to  find  the  consistent  tangent  operator,  the  closest  point  projection  of  the 
stresses  o^,  given  by  Equation  2.19  onto  the  yield  surface  has  to  be  found.  The  yield 
surface  for  the  elliptic  cap  is  characterized  by 


By  taking  derivatives  of  fj  with  respect  to  invariants  and  substituting  into  Equation  2.22 
results 


/2  = (7,-T)2  - (X-L)^  + 


(233) 


I 


= - 18t:aa<<1)(^2)) 


(2.34a) 
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/y. 


2./I+1 


1 ^IGLXiWi)) 


(2.34b) 


Here  we  have  four  unknowns;  I,  / J2„^„  ((|)(f2))  and  L„^,.  Equation  2.13  relates  to 

((J)(f2)),  Equation  2.12  gives  L„^,  in  terms  of  X„^,  and  (<l)(f2))  is  given  by  Equation  2.5.  So 
in  total  there  are  five  equations  with  five  unknowns  and  the  other  three  equations  look 
like 


n*\ 


D 


Ln 


1- 


^vn 


W 


(234c) 


■^n*l 


- 


Rye^'^""  - RdL 


n*l 


(234d) 


(W2))  = 


/o 


N 


(2.34e) 


The  equations  can  not  be  combined  into  one  explicit  equation  for  one  unknown 
because  of  their  nonlinear  form.  In  order  to  solve  the  five  unknowns  we  have  an  inner 
Newton  iteration  scheme  (for  solving  L„+| ) and  an  outer  Newton  iteration  scheme  to 
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solve  I In  order  to  do  this  all  variables  are  expressed  in  terms  of  L„+,  and  I as 
follows.  From  Equation  2.34a,  we  get 


18/: 


(2.35) 


Substituting  Equation  2.35  into  X,,^, , results  in 


X. 


/i+i 


D 


Ln 


3/:^T-3/:e^ J 


7,KW 


Plugging  X„^i  into  Equation  2.34  d,  results  in 


— Ln 
D 


3xif-3/:6:-(/:,-/„^,) 


^KW 


(236) 


again  using  Equation  2.35,  Equation  2.34b  results  in 


/y. 


2./I  + 1 


I ^ n*\  w*!-' 

9X(7  ,-L  ,) 

' /i  + l n*\^ 


(237) 
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The  inner  Newton  iteration  for  L„+,  is  obtained  from  Equation  2.36, 

rk*\  _ j k 

^n*\  ~ ^n*\  ~ ; (2.38a) 


and  G '(!„*,)  = 1 


The  outer  Newton  iteration  for  I„+,  is  obtained  from  Equation  2.34a, 


*+l  jk  ^(4i  + l) 


/ = I - 


ff'iCo 


(238b) 


where 


sCi 


18/fAA.a(4)(^)) 


1- 


dL 


n*\ 

dCi 


aL^i  /ai„+,  and  a(4>(f2))  /ai^,  are  obtained  from  Equations  2.36  and  2.34e  respectively,  in 
the  process  we  need  partial  derivatives  of  fj  with  respect  to  I„^|,  / J2  „^, , and  L^,  which 
can  be  obtained  from  Equation  2.33.  d/ J2_^,  /ai„^,  is  obtained  from  Equation  2.37. 

Once  the  closest  point  projection  of  the  trial  stress  onto  the  yield  surface  o^,  is 
found , the  next  thing  is  the  determination  of  consistent  tangent  operator.  The  tensorial 
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form  of  the  operator  is  the  same  as  that  given  in  Equation  2.28.  The  derivatives  of  dl, 
/5e„^, , and  d^T are  evaluated  from  Equation  2.34a  and  2.34b  and  the  final 
forms  are 


dl. 


n*\ 


de 


n + 1 


3/:.  1>KB  18/i:AA(/^^,-T^^,)£)3  ^ 


i8/:ax(/„.,-z„,,)d, 


a/y. 


2jt*\ 


de 


/i+i 


b: 


1 + 


/2G5,  - 


where 


D2=(1+B  + B/D,) 
B=18KAA.  {({iCQ) 


B,=(  1 + 2 G R'  AA  ((jiCQ) ) 
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Oj  = jiWi)} 


3A:C,  30C, 


3/:C  7>KC. 


D^*Di 


D, 


3OC3 

Df.D, 


D,  = y <<Kf2»  ^ 


c. 


C| , C2 , and  C3  are  partial  derivatives  of  fj  with  respect  to  , and  L^, 

which  can  be  obtained  from  Equation  2.33  . 

Finally  the  consistent  tangent  operator  for  the  elliptic  cap  is 


da 


de 


«+l 


3K  3KB 


D2  Z)|  *^2 


D.-2Gi 


1 ® 1 


rt  ® 1 - 


1 ® n 


b: 


(2.39) 


2G 


B, 


b: 


D,-2Gi 


n®n  + 2G^I 


where5=/J2^,//J2  Vi  • 

Tension  Cutoff:  The  yield  surface  for  tension  cutoff  is  characterized  with  ideal 


plasticity  as 
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f,-T  - and  "9^AA(({)(/^)), 
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The  Newton  algorithm  was  used  to  solve  I^+i , 


jk*\  _ jk  _ ^(4  + i) 
^n*\  ~ ^n*\ 


^(C.) 


(2.41) 


where 
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To  evaluate  the  consistent  tangent  operator,  Equation  2.28,  the  partial  derivates, 
5Ii  n+,  /de^i , and  5/ J2.n+i/^€n+i  niust  be  determined: 
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and  the  final  consistent  tangent  operator  is  given  as 
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(2.43) 


Tension  corner  region  is  shown  in  Figure  1.1  and  is  characterized  with  f,  > 0, 

^ 0,  A(j)  I > 0,  A4>3  > 0 , and  A(j)2  = 0.  Since  both  the  failure  and  the  tension  cutoff 
surfaces  behave  as  in  ideal  plasticity,  the  return  point,  o„,.„  is  at  the  intersection  off,  and 
fj,  i.e., 

/,  „,!  = r , and  - F^(T)  (2.44) 


Simo  et  al.  (1988a)  found  the  inequalities  that  characterized  this  region  based  on  closest 
point  projection  to  be 
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(2.45) 


To  evaluate  the  consistent  tangent  operator.  Equation  2.28,  the  partial  derivates, 

^ must  be  determined.  The  latter  are  evaluated  from  Eqs. 
2.45  along  with  the  consistency  conditions  df|/d€„^,  = 0 , and  dfj/de^,  = 0 . From 
af3/ae„,i  = 0 and  Equation  2.44,  it  can  be  shown  that  dl,  = 0,  and  from  af,/ae„^, 

= 0,  then  a/ J2,„+i/ae„+|  = 0.  The  final  form  of  the  consistent  tangent  operator  is 
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(2.46) 


Compressive  comer  region  is  shown  in  Figure  1 . 1 and  is  characterized  with  f,  > 0, 
> 0,  A(|)|  > 0,  A(|),  > 0 , and  A(J>3  = 0.  From  the  assumed  hardening  law  it  follows 
that  the  failure  surface  and  cap  behave  as  in  ideal  plasticity  in  this  region.  Therefore,  the 
return  point,  o^,,  is  at  the  intersection  of  f,  and  fj,  i.e., 

~ ^e(^n)  (2.47) 


Simo  et  al.  (1988a)  found  the  inequalities  that  characterized  this  region  based  on  closest 
point  projection  to  be 
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(2.48) 


The  consistent  tangent  operator.  Equation  2.28,  with  the  partial  derivates,  51,  „+,/8e^„ 
and  must  be  determined.  From  Equation  2.48  and  the  consistency 

conditions  6f,/3e^,  = 0 , and  df2lde„+^  = 0 , it  can  be  shown  that  5I,^,/5€n+,  = 0,  and 
^ J2,n+/^Sn+i  ~ ® which  results  in  the  same  consistent  tangent  operator  as  given  in 
Equation  2.46 
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2.5  Validation  of  the  Viscoplastic  Finite  Element  Code 

The  viscoplastic  algorithm  implemented  in  the  general  purpose  nonlinear  finite 
element  program  PLASFEM  is  taken  through  several  examples  as  mentioned  in  Katona 
(1984).  The  problems  are  modelled  using  a single  brick  element  (unit  cube)  with 
appropriate  boimdary  conditions.  The  first  in  the  series  is  McCormic  Ranch  Sand 
example,  which  is  used  to  illustrate  the  behavioral  aspects  of  the  viscoplastic  model,  like 
the  effect  of  fluidity  parameter  on  the  stress  response  for  the  same  strain  loading  history. 
To  illustrate  the  capabilities  and  limitations  of  the  model,  sets  of  experimental  data  are 
considered,  such  as  (hard  limestone,  sedimentary  rock,  well-graded  sand).  The  model  is 
also  tested  against  wide  range  of  loading  paths  and  loading  rates. 

In  order  for  the  model  to  perform  on  par  with  the  experimental  data,  the  model 
parameters  have  to  be  identified  properly.  For  detailed  identification  of  parameters  refer 
Katona  (1984). 

2.5.1  McCormick  Ranch  Sand 

In  this  example,  a uniaxial  stress  response  is  assessed  under  uniaxial  strain 
loading.  The  viscoplastic  cap  model  employs  the  plasticity  parameters  for  McCormick 
Ranch  Sand  given  by  Sandler  and  Rubin  (1979).  The  loading  history  and  model 
parameters  are  given  in  Figure  2.1.  The  model  is  assessed  for  different  fluidity 
parameters  to  see  the  effect  on  the  stress  response.  The  time  steps  used  were  At  = 0.01, 
0.1  and  the  results  obtained  were  identical.  Figure  2.2  shows  computed  values  of  axial 
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stress  o„  vs  time  for  different  values  of  the  fluidity  parameters,  i.e.,  y=0.001,0.01,  & 0.1. 
For  large  values  of  y=0.1,  the  response  is  close  to  inviscid  (elasto-plastic)  throughout  the 
loading  history,  this  is  because  the  viscous  strain  occurs  very  rapidly.  It  should  also  be 
observed  that  the  results  corresponding  to  steady  state  values  of  o„  agree  with  the 
inviscid  plasticity  solution,  equation  2.4,  thereby  validating  the  viscoplastic  model.  For 
any  value  of  y>0. 1 , the  solution  resembles  the  inviscid  solution  over  the  entire  range  this 
is  because  the  small  peaks  at  time  = 1.0  and  5.5  will  flatten  further.  When  y is  reduced, 
i.e.,  when  y=0.001,  the  greatest  magnitude  of  peak  stress  occurs,  this  is  because  the 
viscoplastic  strain  rate  is  reduced  and  elastic  strain  rate  is  increased.  The  response  will 
approach  purely  elastic  when  y— ► 0. 

2.5.2  Limestone  in  Triaxial  Stress 

Robertson  (1960)  conducted  an  elaborate,  non  standard,  triaxial  test  experiments 
on  specimens  of  Solenhofen  limestone  specimens  to  measure  axial  strain  €„  resulting 
from  a variable  axial  stress  loading  sequence.  The  loading  history  and  the  parameters  are 
shown  in  Figure  2.3.  As  seen  from  the  Figure  2.3  an  initial  triaxial  stresses  of  (o,j  = 

96.1  ksi,  and  Ojj  = Ojj  = 44.1  ksi)  are  rapidly  imposed.  The  lateral  stresses  are  held 
constant  throughout  the  experiment  and  the  axial  stress  is  intermittently  step  loaded  at 
time  = 7.2,  12.9  & 22.8  ksec.  After  each  step  loading,  o,,  decreases  by  some  amount 
(magnitude  is  known  but  not  the  time  history)  due  to  the  nature  of  the  testing  apparatus. 
The  axial  strains  were  reported  before  and  after  each  loading  step. 


Elastic  Moduli:  Bulk  = 66.7  ksi;  Shear  = 40.0  ksi 

Failure  Surface:  a = 0.25  ksi;  P = 0.67  ksi  ';  y = 0-18  ksi 
Cap  Surface:  R = 2.5;  fo=  0.25  ksi;  Xq  = 189.0  ksi 

Cap  Hardening:  W = 0.066  ; Dg  = 0.67  ksi  ' 

Viscous  Flow  Function:  N = 1;  fg  = 0.25  ksi 


Time  (seconds) 


Figiue  2.1  Uniaxial  Strain  Schedule  for  Sand 
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Time  (seconds) 


Figure  2.2  Axial  Stress  Response  for  Three  Values  of  y 


The  test  specimen  was  modelled  using  a single  8 noded  3D  brick  element.  The 
failure  surface  used  was  a standard  Drucker-Prager  form  and  the  initial  cap  siuface  was 
located  well  into  the  compression  region  by  setting  = 212.0  ksi.  This  setting  was 
necessary  to  provide  a large  elastic  region,  so  that  the  initial  jump  loading  does  not  cause 
excessive  viscoplastic  flow.  The  loads  were  applied  as  follows:  the  initial  axial  load  of 
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96.1  ksi  was  applied  within  100  seconds  in  5 steps  and  the  lateral  load  is  applied  at  the 
first  time  step  in  full  i.e.,  44.1  ksi)  and  at  every  step  the  load  is  assumed  to  be  applied 
gradually  in  3 time  steps  over  a span  of  75  secs.  Figure  2.4  shows  the  strain  history  data 
points  along  with  the  viscoplastic  cap  model  representation.  The  model  produces  a fairly 
good  correlation.  It  is  also  observed  that  no  viscoplastic  flow  has  occurred  during  the 
jumps.  This  was  best  achieved  by  having  a constant  bulk  modulus  and  variable  shear 
modulus  monotonically  decreasing  with  -f  J2 . 

2.5.3  Sedimentary  Rock  in  Triaxial  Stress 

AJcai  et  al.  (1977)  have  conducted  standard  triaxial  tests  to  investigate  the 
viscoplastic  yielding  of  soft  sedimentary  rock.  The  data  presented  here  consists  of  four 
separate  creep  tests  with  the  same  confining  pressure  but  with  different  axial  loads.  For 
each  test,  the  axial  load  is  rapidly  applied  and  held  constant.  The  duration  of  the  test  is 
about  8000  min.  Figure  2.5  defines  imposed  stress  states,  initial  cap  setting,  and 
parameters  for  soft  sedimentary  rock. 

The  elastic  properties  were  determined  such  that  for  each  axial  load  imposed,  the 
initial  strains  were  elastic  and  no  viscoplastic  flow  occurs.  The  failure  surface  used  is  a 
standard  linear  Drucker-Prager  surface  and  the  initial  cap  setting  was  taken  well  into  the 
compression  range  by  setting  = 800  psi.  The  motivation  behind  these  choices  were 
that  test  1 exhibited  pure  elastic  behavior,  test  2 creep  and  steady  state,  and  test  3 & 4 
showed  secondary  creep  and  tertiary  creep. 
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Elastic  Moduli:  Bulk  = 3500  ksi; 

Shear,  G(  J2 ) = 91(l+21exp(-0.0012  Jj  ) ) ksi 
Failure  Surface:  a =1.0  ksi;  0 = 0.25 
Cap  Surface:  R = 2.4;  fo  = 1 .0  ksi;  Xq  = 2 1 2.0  ksi 

Cap  Hardening:  W = 0.55  ; Dq  = 0.0024  ksi  ' 

Viscous  Flow  Function:  N = 1 ; fo  = 1 .0  ksi;  y = 0.2x  10  ’ sec  ' 


Figure  2.3  Triaxial  Stress  Loading  Schedule  and  Model  Parameters  for  Limestone 

(Experimental) 
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Time  x 10^  (seconds) 


Figure  2.4  Axial  Strain  Response  Data  and  Viscoplastic  Model  Representation 
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The  triaxial  test  was  modelled  by  finite  element  program  using  a single  8 nodded 
brick  element.  The  confining  pressure  and  the  axial  load  are  applied  in  a single  time  step 
and  held  constant  throughout  the  analysis.  Four  separate  runs  were  made  for  different 
axial  loads.  Figure  2.6  shows  measured  deviatoric  strain  e„  reported  by  Akai,  et  al., 
superimposed  with  the  finite  element  prediction  by  the  viscoplastic  model. 

The  elastic  response  in  test  1 , creep  & the  steady  state  viscoplastic  response  in 
test  2 are  very  well  matched  with  the  experimental  data.  The  model  is  also  able  to  match 
the  primary  creep  and  secondary  creep  as  observed  in  tests  3 & 4 but  not  the  tertiary  creep 
as  shown  by  the  last  two  points  in  test  3 and  from  the  beginning  in  test  4.  This  can  be 
over  come  by  introducing  strain  softening  into  the  hardening  function. 

2.5.4  Sand  in  Uniaxial  Strain  with  Variable  Load  Rates 

This  ingenious  experimental  test  was  conducted  at  the  Army's  Waterways 
Experimental  Station  (8).  This  test  was  undertaken  to  directly  assess  the  effect  of  loading 
rate  on  the  constitutive  behaviour  of  a dry  remolded  sand  (20-40  Ottawa  Sand).  The  sand 
was  molded  into  a thin  disk-shaped  specimen  at  the  bottom  of  a rigid  cylinderical  test 
chamber  that  provided  lateral  constraint  (uniaxial  strain).  By  means  of  elaborate  ram  and 
explosive  loading  devices,  several  specimens  were  pressure  loaded  with  different  rise 
times  ranging  from  approximately  0.2-20,000  msec.  Data  for  each  test  included  the  time 
histroy  of  the  pressure  loading  and  the  corresponding  strain  history. 


Elastic  Moduli:  K = 125,000  psi;  G = 60,000  psi 
Failure  Surface:  a = 275  psi;  0 = 0. 1 1 
Cap  Surface:  R = 0.35;  Xq  = 800.0  ksi 

Cap  Hardening:  W = 0.25  ; Dg  = 0.00078  psi"' 

Viscous  Flow  Function:  N = 1 .6;  fg  = 275.0  ksi;  y = 0-5x  1 0'^  sec’' 


li  (psi) 


Figure  2.5  Stress  States,  Initial  Cap  Settings  for  Soft  Sedimentary  Rock 


Deviatoric  strain,  e 
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0 10  20  30  t 40  50  60  70  80  90 

Time  x 10^  (minutes) 


Figxire  2.6  Deviatoric  Strain  History  Responses  for  Soft  Sedimentary  Rock 
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Model  Parameters 

Elastic  Moduli:  K = 3530  MPa;  G = 1630  MPa  psi 
Failure  Surface:  a = 5.0  MPa;  P = 0.004  MPa‘‘;  4.5  MPa 
Cap  Surface:  R = 2.0;  X(,  = 3.0MPa 

Cap  Hardening:  W = 0.027  ; Dq  = 0.0 1 4 MPa'' 

Viscous  Flow  Function:  N = 1 .0;  fo  = 1 .0  MPa;  y =2. Ox  10'*  MSec"' 


Figure  2.7  Slow  Loading  Stress  Histories  (Experimental) 
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Figure  2.8  Rapid  Loading  Stress  Histories  (Experimental) 


Figure  2.7  shows  the  pressure  loading  history  (o,,  stress)  for  the  slow  loading  rate 
having  a rise  time  of  15,000  msec.  At  the  other  extreme,  Figure  2.8  shows  the  stress 
histories  for  the  'rapid'  loading  rate  having  a rise  time  of  0.3  msec.  The  problem  was 
modelled  by  finite  element  using  a single  8 nodded  brick  element  with  suitable  boundary 
conditions  to  simulate  experimental  conditions  (uniaxial  strain).  The  resulting  stress- 
strain  curves  for  both  loading  rates  are  shown  in  Figure  2.9.  Here  it  is  clearly  evident  that 
the  sand  specimen  exhibits  rate  dependent  stress-strain  behavior. 


46 


Intermediate  loading  rates  with  rise  times  on  the  order  of  100  msec  were  also 
carried  out  and  the  results  yielded  were  almost  identical  to  the  slow  loading  rate 
experiment.  From  these  findings  two  important  observations  can  be  made:  (l)The 
nonlinear  stress-strain  relationship  for  the  slow  loading  rate  is  not  time  dependent  and 
may  be  assumed  invisicid;  and  (2)  rate  effects  only  become  significant  when  the  rise  time 
approaches  the  sub-millisecond  range.  These  observations  were  used  in  identifing  the 
model  parameters. 

From  Figure  2.9,  it  can  be  concluded  that  the  viscoplastic  cap  model  accurately 
reflects  the  test  data.  For  the  slow  loading  test,  the  cap  surface  continually  moves  with 
the  stress  state  producing  a stiffening  stress-strain  response  as  shown.  Since  the  stress 
state  is  on  the  cap  surface  during  loading,  the  immediate  unloading  response  is  initially 
elastic  prior  to  re-entering  the  failure  surface.  On  the  other  hand,  for  the  fast  loading  rate, 
the  cap  surface  lags  behind  the  loading  stress  state  producing  an  'apparent'  softening 
stress-strain  response  (which  is  really  a time-dependent  effect).  Just  prior  to  unloading, 
the  stress  state  is  well  above  the  cap  surface.  Thus,  when  unloading  occurs,  the  stress 
state  remains  in  the  viscoplastic  domain  producing  additional  strain  accumulations  as 
shown.  The  correlation  between  the  model's  performance  and  the  observed  performance 
is  quite  remarkable,  particularly  \vith  regard  to  matching  the  rapid  loading  behavior, 
which  is  characteristic  of  ground  shock  problems. 
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Figure  2.9  Stress/Strain  Response  for  Slow  and  Rapid  Loading 


CHAPTER  3 

PULLOUT  SIMULATION  OF  POSTINSTALLED  ADHESIVE-BONDED  ANCHORS 

3.1  Introduction 

All  concrete  anchors  in  use  today  may  be  divided  into  one  of  two  categories:  cast- 
in-place  and  postinstalled.  Postinstalled  anchors  include  expansion  anchors,  undercut 
anchors  and  more  recently  adhesive  bonded  anchors.  In  the  past,  the  majority  of  all 
anchors  were  cast-in-place  or  post-installed  expansion  and  undercut  anchors;  however 
with  the  advent  of  high  strength  polyester,  vinylester  and  epoxy  adhesives,  the  use  of 
adhesive-bonded  anchors  has  increased  significantly  over  the  past  few  years.  Typically 
used  on  structural  repairs  or  retrofit,  adhesive-bonded  anchors  are  now  being  used  on  new 
structures  due  to  diminished  cost,  time,  and  ease  of  installation. 

A adhesive-bonded  anchor  is  a reinforcing  bar  or  threaded  rod  inserted  into  a 
drilled  hole  (usually  10-25%  larger  than  the  diameter  of  the  anchor)  within  hardened 
concrete  with  a structural  adhesive  acting  as  a bonding  agent  between  the  concrete  and 
steel  anchor.  Current  design  often  employs  the  use  of  proprietary  design  tables  provided 
by  manufacturers.  These  tables  are  usually  based  on  average  strengths  determined  from 
tests  on  specific  diameters  and  embedment  lengths  of  rebar  or  threaded  rods.  Measured 
strengths  are  then  reduced  by  a factor  to  establish  allowable  anchor  strengths. 
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Eligenhausen  et  al.  (1984)  proposed  an  equation  from  numerous  lab  tests  to 
predict  pullout  force  on  anchors  exhibiting  concrete  cone  failures.  Subsequently,  Cook  et 
al.  (1993)  developed  an  equation  to  predict  bond  failure  for  pullout  tests  from  equilibrium 
and  kinematics  between  the  threaded  rod  and  the  adhesive  layer.  Cook  (1993)  combined 
the  work  of  Eligenhausen  et  al.  (1984)  and  Cook  et  al.  (1993)  to  predict  the  pullout  force 
and  cone  depths  of  anchors  in  the  100  mm  to  200  mm  range.  Both  of  the  latter 
researchers'  expressions  for  pullout  resistance  require  an  estimate  of  either  the  average  or 
maximum  shear  stress  within  the  adhesive  bond  layer. 

Finite  element  studies  of  various  proprietary  anchor  designs  have  been  reported  by 
Peier  (1983)  and  by  James  et  al.  (1987).  Peier  studies  were  limited  to  short  anchors  with 
embedment  length  to  diameter  ratios  of  3.2  . The  James  et  al.  (1987)  study  covered  a 
number  of  embedment  lengths  to  diameter  ratios,  as  well  as  different  epoxy  properties. 
However  both  studies  involved  linear  elastic  (ANSYS  code)  and  elastic-perfectly  plastic 
(ABAQUS  code)  analyses  of  bonded  anchors  with  bolt  heads  at  the  embedded  end  of  the 
anchors.  Typically,  adhesive-bonded  anchors  do  not  have  a bolt  head  or  nut  at  the 
embedded  end  of  the  anchor. 

This  chapter  reports  on  a state  of  the  art  elasto-plastic  finite  element  simulation  of 
pullout  of  postinstalled  anchors  with  embedment  depths  to  anchor  diameter  ratios 
varying  from  4 to  40.  The  simulation  revealed  that  the  concrete  failure  cone  was  found  to 
initiate  at  the  concrete-adhesive  bond  interface  and  progress  towards  the  surface  as  a 
tension  zone.  Experimental  testing  of  postinstalled  anchors  validated  this  concept  as  well 


as  substantiating  the  predicted  pullout  forces  for  a number  of  embedment  to  anchor 
diameter  ratios. 
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3.2  Background 

Adhesive-bonded  anchors  transfer  the  load  from  the  steel  rod  through  the  adhesive 
layer  into  the  concrete  along  its  bonded  interface.  Figure  3.1  displays  four  types  of  failure 
modes  for  bonded  anchors  (Cook,  1993): 

(a)  Anchor  steel  failure:  This  is  characterized  by  yielding  and  subsequent  fracture 
of  the  steel  rod;  usually  occurs  for  long  embedment  lengths,  Figure  3.1  (a); 

(b)  Concrete  cone  failure  only:  As  noted  in  Luke  et  al.  (1985),  this  is  observed  in 
only  shallow  embedments  (75  mm  or  smaller).  Figure  3.1  (b); 

(c)  Bond  failure:  This  may  occur  if  the  bonded  surface  lacks  adequate  strength. 
Cook  et  al.  (1993),  shows  that  this  failure  can  occur  when  the  top  portion  of  the 
embedment  length  is  debonded  about  50  mm.  Cook  et  al.  (1991)  produced  this  failure  by 
performing  confined  tension  tests.  Figure  3.1  (c); 

(d)  Combined  cone-bond  failure:  For  embedments  greater  than  50-100  mm,  the 
commonly  observed  failure  is  by  the  combined  cone-bond  failure  mode,  with  a shallow 
cone  attached  to  the  top  of  the  anchor.  Figure  3.1  (d). 

Luke  et  al.  (1985)  reported  that  bond  failure  occurred  first  followed  by  cone 
failure.  Subsequently  Cannon  et  al.  (1981)  reported  that  the  shallow  cone  formed  prior  to 
failure  of  the  anchor.  Since  then  Collins  et  al.  (1989)  measured  displacements  at  the 
surface  of  the  concrete  adjacent  to  anchor  and  at  embedded  end  of  anchor  and  noted  that 
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Figure  3.1  Typical  modes:  (a)  Steel  Failure;  (b)  Concrete  Cone  Failure;  (c)  Bond 
Failure;  (d)  Combined  Cone-Bond  Failure 
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the  cone  and  bond  failure  occurred  simultaneously.  The  Cook  et  al.  (1993)  results  also 
report  a combined  cone  and  bond  failure  occurring  together,  as  well  as  stating  that  the 
load-deformation  response  of  the  anchor  was  imreliable  once  bond  failure  occurred.  The 
erratic  behavior  after  bond  failure  was  postulated  to  be  due  to  mechanical  interlock 
between  adhesive  layer  and  concrete  interface. 

Currently,  there  is  no  accepted  way  to  accurately  determine  the  strength  of 
adhesive-bonded  anchors.  Several  behavioral  models  have  been  proposed  to  determine 
the  tensile  strength  of  adhesive-bonded  anchors  based  on  their  failure  modes  (Cook  et  al., 
1993,  Cook,  1993).  A number  of  methods  will  be  reviewed  for  each  of  the  identified 
failure  methods. 

3.2.1  Concrete  Cone  Model 

Based  on  tests  of  adhesive-bonded  threaded  rod  anchors,  Eligehausen  et  al. 
(1984),  suggested  the  following  equation  for  pullout  capacity,  where  N„  is  the  failure 


pullout  force,  h^f  is  the  embedded  depth  of  the  anchor,  and  f^  is  the  compressive  strength 
of  the  surrounding  concrete  at  time  of  testing.  Rehm  et  al.  (1988)  compared  Equation  3.1 
with  test  results  for  polyester  bonded  anchors  with  different  embedment  lengths  and 
diameters  (d„ ) with  the  same  h,/do  ratios.  The  correlation  between  predicted  and 


{newtons) 


(3.1) 
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measured  was  quite  good;  however  the  method  may  be  inappropriate  for  h^/do  ratios  other 
than  the  one  tested. 

3.2.2  Bond  Models 


different  shear  stress  profiles  within  the  adhesive  layer.  For  the  uniform  bond-stress 
model,  the  pullout  resistance,  N„  is  simply  given  as 


where  t,  is  the  uniform  shear  stress  in  the  adhesive  layer,  d„  is  the  diameter  of  the  hole, 
and  hjf  is  the  embedded  anchor  depth.  For  the  elastic  bond  stress  model  proposed  by 
Cook  et  al.  (1993),  the  pullout  resistance,  N„  is  given  as 


where  X’  is  an  experimentally  determined  value  for  stiffness  of  the  anchor  that  includes 
the  effects  of  the  elastic  shear  modulus  of  the  adhesive  layer  and  the  axial  stififiiess  of  the 
threaded  rod,  d„  is  the  diameter  of  the  hole,  and  is  the  maximum  shear  stress  in  the 
adhesive  layer.  It  was  found  that  the  uniform  stress  model  and  the  elastic  bond-stress 
model  give  similar  results  for  embedment  depths  up  to  40/ do  (Cook,  1993)  after  which 
the  uniform  model  results  in  an  over  prediction  of  the  pullout  resistance. 


There  are  two  bond  models  which  have  been  proposed  for  design.  Both  employ 


(3.2) 


(3.3) 
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3.2.3  Combined  Cone-bond  Model 


For  postinstalled  anchors  with  embedment  to  diameter  ratios  greater  than  5, 
combined  cone  and  bond  failure  are  known  to  occur.  Cook  (1991)  combined  the  pullout 
resistance  of  the  cone  and  bond: 


N =N  +7V,  , 

u cone  bona 


(3.4) 


and  minimized  it  with  respect  to  cone  depth.  Using  the  Eligehausen  et  al.  (1984) 
expression  for  cone  failure  with  the  uniform  bond-stress  model  and  minimizing  with 
respect  to  h^„j , Cook  (1993)  arrived  at  pullout  resistance  of 

= 0.92h^onc)lfc^'^o'^dXh^f-h^^)  (newtons)  (3.5) 


In  the  minimization,  the  cone  depth,  h,^  , was  found  to  be 

h ~ 

^*conc 

1.84 


(3.6) 


For  the  elastic  bond-stress  model  of  Cook  et  al.  (1993)  and  the  Eligehausen  et  al.  (1984) 
expression  for  cone  failure  (Cook,  1993)  arrived  at  a ultimate  pullout  resistance  of 


(3.7) 


and  a failure  concrete  cone  depth  of 
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{mm) 


(3.8) 


The  only  difference  between  Equations  3.6  and  3.8  is  the  value  of  Tq  vs  and 
the  inclusion  of  the  bracketed  term  in  Equation  3.8.  For  short  embedments,  and 
are  very  similar  (Cook,  1993)  and  the  bracketed  term  in  Equation  3.8  is  close  to  unity.  In 
the  case  of  deeper  embedment  depths  is  greater  than  and  the  bracket  term  in 
Equation  3.8  is  less  than  unity  resulting  in  smaller  failure  cone  depths  from  Equation  3.8. 
Experimental  testing  by  Collins  et  al.  (1989)  and  Doerr  et  al.  (1989)  revealed  a similar 
trend.  Both  Equations  3.6  and  3.8  require  the  user  to  provide  and  from  either 
partially  bonded  or  confined  anchor  pullout  tests  in  concrete  of  similar  strengths  and 
adhesive  bonding  agents. 

3.2.4  Elasto-Plastic  Finite  Element  Model 

The  model  used  herein  for  both  the  concrete  and  epoxy  adhesive  was  the 
Sandler  et  al.  (1976)  cap  model.  Simo  et  al.  (1988a)  showed  the  Sandler  et  al.  (1976) 
model  did  an  excellent  job  of  predicting  concrete  response  for  a multitude  of  different 
stress  paths.  For  constitutive  model  details  refer  chapters  1 and  2.  A complete 
description  of  the  material  parameters  used  to  model  the  concrete  and  adhesive  bond  layer 


follows. 
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3.3  Experimental  Program 

The  following  describes  the  experimental  testing  that  was  performed  to  validate 
the  analytical  model.  The  anchors  were  postinstalled  in  0.38  m deep  by  1.22  m wide  and 
2.4m  long  concrete  slabs.  The  slabs  were  cast  from  ready-mix  Florida  Department  of 
Transportation  (FDOT)  Class  II  concrete  with  a minimum  specified  design  compressive 
strength  of  24.8  MPa  (3600  psi)  at  28  days.  Test  cylinders  were  made  and  cured  under 
the  same  atmospheric  conditions  as  the  slabs.  The  cylinders  were  subsequently  tested 
after  90  days  when  the  installed  anchors  were  to  be  loaded.  Compressive  strengths 
ranged  from  39  MPa  (5700  psi)  to  43.4  MPa  (6300  psi)  and  tensile  strengths  varied  from 
2.5  MPa  (375  psi)  to  3.3  MPa  (475  psi). 

The  anchors  tested  were  ASTM  A 193  Grade  B7  threaded  rods  with  a 15.9  mm 
(5/8  in.)  diameter.  The  minimum  tensile  strength  of  the  rods  were  862  MPa  (125  ksi)  to 
ensure  that  adhesive  bond  failure  would  occur  before  the  steel  began  to  yield  under 
loading.  The  embedment  depths  of  the  anchors  varied  from  76  mm  to  152.4  mm. 

The  anchor  holes  were  drilled  with  a rotary  hammer  with  a 19  mm  (3/4  in.) 
diameter  drill  bit.  The  holes  were  cleaned  using  a stiff  bristle  brush  and  a compressed  air. 
All  anchors  were  installed  in  dry  holes  as  recommended  by  the  manufacturer.  Care  was 
exercised  to  ensure  that  all  anchors  were  vertically  aligned  and  that  0.355  m of  clear 
distance  existed  between  adjacent  bolts  (i.e.  no  disturbance  from  adjacent  anchors).  The 
adhesive  bond  was  allowed  to  cure  at  room  temperature  for  7 days  before  testing. 
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The  unconfmed  tests  were  performed  in  accordance  with  the  ASTM  E488  (1984) 
specification.  Figure  3.2  shows  the  frame  that  was  designed  and  constructed  to  conduct 
the  tests. 

A total  of  eighteen  fully  bonded  15.9  mm  diameter  anchors  were  tested.  The 
embedment  depths  were  76  mm,  102  mm  127  mm  and  152  mm.  All  of  the  depths  with 
the  exception  of  the  1 52  mm  embedment  were  repeated  5 times  each;  for  the  1 52  mm 
only  3 tests  performed. 

Presented  in  Figure  3.3  is  a picture  of  some  of  the  failed  anchors  with  their 
concrete  cone  attached.  Evident  from  the  Figure  3.3,  the  76.2  mm  embedded  anchors  had 
complete  concrete  cone  failure.  For  all  anchors  with  embedment  depths  greater  than  102 
mm,  a combination  of  concrete  failure  (cone)  and  adhesive  bond-concrete  interface 
failure  occurred.  Note  that,  the  concrete  cones  on  the  failed  anchors  were  well  attached  to 
their  respective  anchors.  Also,  the  size  of  the  concrete  failure  cones  varied  for  a given 
anchor  embedment  depth;  i.e.  the  102  mm  embedded  anchors  had  failure  cone  depths 
ranging  from  41.6  mm  to  71.6  mm,  and  the  127  mm  embedded  anchors  had  failure  cone 
depths  varying  fi-om  37  mm  to  57  mm.  Shown  in  Figure  3.4  is  a concrete  test  slab  with 
postinstalled  anchors,  as  well  as  concrete  failure  cones  from  previous  anchor  tests.  Table 
3.1  shows  the  mean  pullout  resistance  and  range  of  failure  cone  depths  versus  embedment 
depth  for  all  the  anchors  tested.  The  same  epoxy  (amine)  adhesive  was  used  to  bond  all 


anchors. 
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Figure  3.2.  Unconfined  Anchor  Pullout  Testing  Apparatus 
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Figure  3.3  Postinstalled  Failed  Anchors:  76.2  mm,  101.6  mm,  127  mm  and  152.4  mm 

Embedments 


Figure  3.4  Postinstalled  Anchors  and  Failure  Cones  in  Concrete  Test  Slabs 
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Table  3.1  Anchor  Pullout  Test  Results 


Embedment  Depth 

Average  Pullout  Resistance 

Average  Failure  cone  Depth 

76.2  mm 

50  kN 

76.2  mm 

101.6  mm 

72  kN 

40  mm  - 72  mm 

127  mm 

84  kN 

37  mm  - 56  mm 

152.4  mm 

105  kN 

41  mm  - 76  mm 

3.4  Numerical  Analysis 

Pullout  of  adhesive-bonded  anchors  were  simulated  using  the  code  PLASFEM,  an 
elasto-plastic  finite  element  code  under  development  at  the  University  of  Florida.  The 
finite  element  mesh  used  to  model  the  unconfmed  fully  bonded  single  anchors  is  given  in 
Figure  3.5.  Eight  nodded  isoparametric  axisymmetric  elements  were  used  to  represent 
all  materials:  steel  anchor,  epoxy  adhesive  and  surrounding  concrete.  The  complete  mesh 
took  advantage  of  symmetry  and  was  489.6  mm  wide  by  508  mm  deep  (mesh  exhibited 
little,  if  any  boundary  effects).  The  boundary  conditions  were  restraint  of  translations,  in 
the  radial  direction  on  the  left  boundary,  in  the  vertical  direction  on  the  bottom  boundary, 
and  along  radial  & vertical  directions  on  the  right  boundary.  The  mesh  had  324  elements: 
1 8 steel  elements  of  which  one  element  was  extended  above  concrete  surface  to  simulate 
the  experimental  conditions;  1 8 epoxy  elements;  and  288  concrete  elements.  The  aspect 
ratio  of  the  elements  in  the  neighborhood  of  the  epoxy  varied  between  3 to  7 depending 
on  the  anchor's  embedment  lengths.  There  were  a total  of  1850  degrees  of  freedom,  with 
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Figiire  3.5  Finite  Element  Mesh  for  Anchor  Bolt  Pullout  Simulation 
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each  element  employing  3x3  Gauss  quadrature  to  determine  stresses.  The  incremental 
stress  applied  to  the  top  of  the  steel  anchor  was  6.9  MPa,  except  for  the  76.2  mm 
embedded  anchor  case  where  it  was  3.5  MPa.  Each  simulation  (anchor  test)  employed 
from  50  to  100  load  steps,  which  took  1 to  3 iterations  to  converge.  A typical  run 
required  1 0 to  25  minutes  of  CPU  time  on  a 90  Mz  Pentium  PC. 

The  constitutive  model  used  to  simulate  both  concrete  and  epoxy  is  Sandler  et  al. 
(1976)  model  described  earlier.  The  steel  anchor  was  modeled  as  a perfectly  elastic 
material.  Based  on  compression  and  split  tension  data,  the  following  failure  strength 
properties  were  assigned  to  the  concrete,  a=26.61  MPa,y=8.00  MPa,  P=0.064  MPa  ',  and 
0=0.05  (see  Figure  1.1);  i.e.  the  unconfmed  compressive  strength  of  concrete  was  39  Mpa 
and  the  tensile  strength  was  2.93  MPa.  The  initial  size  of  the  elliptic  hardening  cap  for 
the  concrete  was  selected  as  L=1 10.32  MPa,  with  an  aspect  ratio,  R=4.43  and  plastic 
hardening  parameters,  W=0.42  and  0=2.21x10  * MPa  from  Simo  et  al.  (1988a) 
recommended  values.  The  elastic  Young's  modulus  of  the  concrete  was  2.83x10'*  MPa 
with  a Poisson's  ratio  of  0. 18.  Simo  et  al.  (1988a)  published  this  data  based  on  well 
documented  results  at  the  University  of  Colorado.  The  program  consisted  of  six  major 
series  of  nonconventional  multi  axial  stress-strain  curves.  The  total  number  of 
experiments  were  67.  One  test  out  of  each  of  the  six  major  series  were  used  to  fit  the 
data.  A total  of  six  tests  were  used  in  the  actual  fit  of  the  model  parameters.  After  the 
optimal  model  parameters  are  obtained,  the  resulting  cap  model  is  used  to  predict  the 
response  of  each  Colorado  test  not  included  in  the  optimization  process  (total  number= 
61).  Considering  the  experimental  data  scatter,  the  predicted  response  is  in  good 
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agreement  with  the  experimental  results.  Since  the  strength  of  the  concrete  used  at  the 
University  of  Florida  was  in  the  same  range  (6300  psi),  the  above  mentioned  parameters 
did  an  excellent  job  in  predicting  the  pullout  capacity  of  an  adhesive-bonded  anchors. 

The  epoxy's  strength  was  characterized  with  a=8.96  Mpa,  y=0.0,  p=0.0,  and  0=0.02  (see 
Figure  1.1);  the  epoxy's  tensile  strength  was  selected  as  10.34  MPa . These  values  were 
obtained  from  the  manufacturer.  Since  there  was  no  published  hardening  data  on  the 
epoxy,  it's  elliptic  cap  was  not  employed  in  the  simulations.  The  steel  anchor  bolt  was 
modeled  with  an  elastic  Young's  Modulus  of  0.2x10®  MPa  and  a Poisson's  ratio  of  0.30. 

In  order  to  see  which  parameters  effect  the  capacity  of  adhesive  bonded  anchors,  a 
parametric  study  was  made  on  the  following  parameters:  different  embedment  depths, 
different  anchor  diameters,  different  glue  thicknesses,  different  concrete  strengths  and 
different  adhesive  strengths. 


3.5  Discussion  of  Results 
3.5.1  Load-Displacement  Characteristics 

Figures  3.6-3. 8 shows  typical  load-displacement  diagrams  obtained  from 
numerical  simulation  for  unconfined  tension  tests  of  bonded  anchors.  The  displacement 
is  measured  at  the  top  of  the  anchor.  The  embedment  lengths  range  from  5 1 mm  to  508 
mm  for  diameters  12  mm,  16  mm  and  25  mm.  The  load  displacement  diagrams  indicate 
the  following:  (a)  the  initial  slope  of  the  load  displacement  diagrams  is  the  same  for  a 
particular  diameter;  (b)  well  defined  inelastic  load  displacement  response  for  all 
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embedment  depths  except  508  mm;  © the  ultimate  load  increases  proportionally  with  an 
increase  in  embedment  depth  up  to  356  mm.  Figures  3.9  and  3.10  support  statement  (c), 
which  show  a plot  of  ultimate  load  vs  various  embedment  depths  for  glue  thicknesses  of 
1.5  mm,  12  mm  and  25  mm.  It  can  be  seen  from  figures  3.9  and  3.10  that  the  increase  in 
load  is  no  longer  linear  at  508  mm.  Figure  3.11  shows  that  the  load  is  fairly  linear  even  at 
508  mm,  which  is  a h^f  /d  of  20.  Since  the  analyses  were  not  run  for  lengths  between  356 
mm  and  508  mm  for  any  of  the  anchor  diameters  under  study,  it  cannot  be  concluded  that 
proportionality  changes  at  356  mm  for  any  diameter.  On  the  other  hand  if  the  results 
were  compared  in  terms  of  h^f  /d,  it  can  be  found  that  the  load  increases  proportionally  up 
to  the  range  of  25  to  30  and  than  drops.  This  can  be  attributed  to  the  anchor  failing  before 
the  adhesive  has  mobilized  its  full  shear  capacity.  Figure  3.12  also  shows  that  the  shear 
stress  mobilized  at  the  interface  is  nearly  constant  at  about  13  MPa  from  a h^f /d  ratio  of 
2 to  about  hjf  /d  of  25  and  then  drops.  This  clearly  shows  that  the  anchor  has  failed 
before  mobilizing  the  full  shear  capacity.  Similar  behaviors  were  observed  for  anchor 
diameters  of  12  mm  and  25  mm. 

3.5.2  Effect  of  Glue  Thickness  on  Ultimate  Load 

The  analysis  was  carried  out  for  glue  thicknesses  ranging  from  1 .5  mm  to  25  mm. 
Figure  3.13  shows  a plot  of  ultimate  load  vs  d^/d  for  an  anchor  of  16  mm  diameter.  The 
change  in  ultimate  load  is  marginal,  less  than  10%  for  a specific  embedment  length  for 
different  glue  thicknesses.  This  trend  was  observed  in  all  embedment  lengths  for 
different  anchor  diameters  such  as  12  mm  and  25  mm.  This  behavior  leads  to  the 
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Figure  3.6  Load  Displacement  Diagram  for  12  mm  Diameter  Anchor  Bolt  for  Various 
Embedment  Lengths  (5 1 mm  - 508  mm) 
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Figure  3.7  Load  Displacement  Diagram  for  16  mm  Diameter  Anchor  Bolt  for 
Various  Embedment  Lengths  (51  mm-508  mm) 
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Figure  3.8  Load  Displacement  Diagram  for  25  mm  Diameter  Anchor  Bolt  for  Various 

Embedment  Lengths  (5 1 mm-508  mm) 
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Figure  3.9  Ultimate  Load  vs  Embedment  Length  Diagram  for  12  mm  Diameter  Anchor 

Bolt  for  Various  Glue  Thicknesses 


Ultimate  Load  ^(,  (kN) 


69 


350 

300 

250 

200 

150 

100 

50 

0 


t = 25  mm 


0 100  200  300  400  500  600 


Embedment  Length  (mm) 


Figure  3.10  Ultimate  Load  vs  Embedment  Length  Diagram  for  16  mm  Diameter  Anchor 

Bolt  for  Various  Glue  Thicknesses 
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Figure  3.11  Ultimate  Load  vs  Embedment  Length  Diagram  for  25  mm  Diameter  Anchor 

Bolt  for  Various  Glue  Thicknesses 
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Figure  3.12  Variation  of  Shear  Stress  for  Different  Embedment  Lengths  for  16  mm 
Diameter  Anchor  Bolt  for  Glue  Thicknesses  of  1.5  mm  and  25  mm 
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Figure  3.13  Load  Vs  Glue  Thicknesses  for  Different  Embedment  Lengths  (51  mm- 
508  mm)  for  16  mm  Diameter  Anchor  Bolt 


conclusion  that  the  shear  capacity  of  the  adhesive  plays  an  imp>ortant  role.  Once  the  shear 
capacity  of  the  adhesive  is  reached  , the  anchor  fails  in  bond  failure  at  the  anchor  / 
adhesive  interface. 

3.5.3  Effect  of  Concrete  Strengths  and  Adhesive  shear  Strength 

A detailed  study  was  done  on  different  concrete  strengths  (20  Mpa,  41  MPa  and  62 
Mpa).  The  analysis  was  repeated  for  different  embedment  depths  ranging  from  51  mm  - 
279  mm  and  different  anchor  diameters  of  12  mm,  16  mm  and  25  mm.  The  failure  load 
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remained  more  or  less  the  same  which  can  be  seen  from  Figures  3.14.  Additional 
calculations  were  performed  with  reduced  shear  capacity  of  the  adhesive  and  it  was  found 
that  the  capacity  of  the  anchor  was  reduced  in  proportion  to  the  reduction  in  the  adhesive 
shear  strength,  Table  3.2.  This  leads  to  the  conclusion  that  the  failures  were  due  to  the 
shear  failures  of  the  adhesive  rather  than  concrete.  These  statements  back  the  findings  in 
the  previous  section. 

At  the  University  of  Florida,  20  different  products  were  tested  by  Cook  et  al. 
(1994)  in  confined  pullout  tests  in  three  different  strength  of  concrete.  The  range  of 
tested  concrete  was  from  27  MPa  to  62  MPa  compressive  strengths.  The  results  of  these 
tests  was  that  the  influence  of  concrete  strength  on  the  anchor  capacity  is  negligible. 


Table  3.2  Shear  Stress  vs  Ultimate  load  for  16  mm  Diameter  Anchor  Bolt 


^eff 

Ultimate  Load  kN 

(mm) 

T = 4.5  MPa 

T = 9.0  MPa 

76 

52 

27 

102 

64 

32 

127 

82 

41.5 

152 

104 

52 

178 

111 

55 
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Figure  3.14  Influence  of  Concrete  Strength 


3.5.4  Failure  Planes  and  Cone  Depths 

Failure  planes  and  depths  to  failure  initiation  were  determined  for  16  mm 
diameter  anchor  with  embedment  lengths  of  76, 102, 127  and  152  mm. 

76  mm  embedment  anchors  The  finite  element  modeling  (Figure  3.15)  revealed 
that  the  76  mm  embedded  anchors  experienced  full  cone  failure.  The  points  in  the  figure 
indicate  the  gauss  points  which  have  reached  either  tension  or  shear  limiting  values.  The 
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points  in  the  shaded  region  have  reached  limiting  value  in  tension  whereas  the  points 
alone  (below  shaded  area)  have  reached  limiting  value  in  shear.  The  hatched  area  is  the 
variation  in  the  experimentally  observed  failure  planes  with  their  corresponding  radii  (124 
mm  and  1 64  mm).  The  distances  h^f  and  h in  Figure  3.15  represent  the  embedment  depth 
of  the  anchor  and  the  depth  of  the  concrete  cone  failure,  respectively.  Note  the  growth  of 
the  tension  limiting  zone  initiated  at  the  epoxy-concrete  interface  at  the  bottom  of  the 
displayed  tension  region  and  progressed  toward  the  surface.  The  dark  circle  (73  mm  from 
the  top  of  the  mesh ) on  Figure  3.15  represents  the  gauss  point  at  which  the  element  stress 
did  not  converge.  Since  the  stresses  depend  on  the  volumetric  strain,  this  can  be  assumed 
as  the  point  where  the  failure  initiated  and  followed  the  line  mentioned  as  numerical. 
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Figure  3.15  Development  of  Failure  in  the  Vicinity  of  the  76  mm  Embedded  Anchor 
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102  mm  embedment  anchors  Displayed  in  Figure  3.16  and  3.17  is  the 
experimentally  observed  concrete  cone  failure  for  the  102  mm  embedded  anchors.  Figure 


3.16  shows  the  plan  area  and  Figure  3.17  displays  the  cross-sectional  view  of  the 
failures.  From  Figure  3.17  it  is  quite  evident  that  the  depth  of  concrete  failure  cone 
varies,  but  the  variation  is  minimal.  It  can  also  seen  from  Figure  3.17  that  the  failure 
cone  consisted  of  two  parts:  a primary  cone  which  originates  from  the  anchor  and  a 
secondary  cone  in  the  vicinity  of  the  concrete  surface.  It  is  postulated  that  the  secondary 
cone  formation  takes  place  when  rigid  body  movement  of  the  adhesive  the  zone  of 
tension  failure,  and  the  dots  represent  shear  failure  of  the  concrete.  Again  the  tension 
zone  originated  at  the  adhesive-concrete  interface  in  the  bottom  of  the  shaded  area  and 
propagates  the  surface.  The  gauss  point  at  which  the  element  stress  did  not  converge  was 
about  62  mm  from  the  top  and  displayed  a 56°  from  the  horizontal  is  the  proposed  failure 
plane  from  the  numerical  analysis.  Also  displayed  in  Figure  3.18  is  the  measured 
variation  in  concrete  failure  cone  depths  from  experimental.  The  analytical  prediction 
was  in  the  middle  of  measured  failure  planes  range. 

127  mm  embedment  anchors  Presented  in  Figure  3.19  and  3.20  is  the 
experimentally  observed  concrete  cone  failure  for  the  127  mm  embedded  anchors.  Figure 
3.19  shows  the  plan  area  and  Figure  3.20  displays  the  cross-sectional  view  of  the  failures. 
The  variation  the  failure  cone  depths  were  37  mm  to  56  mm,  the  primary  cone  radii  were 
58  mm  to  130  mm  and  the  inclination  of  the  primary  cone  failures  were  33°  to  24°  from 
the  horizontal,  respectively.  Displayed  in  Figure  3.21  is  the  finite  element  simulation  of 
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Figure  3.16  Plan  View  of  Concrete  Cone  Failure  for  102  mm  Embedded  Anchor 


Figure  3.17  Cross-sectional  View  of  Concrete  Cone  Failure  for  1 02  mm  Embedded 

Anchor 
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r=90mm 


Figure  3.18  Development  of  Failure  in  the  Vicinity  of  the  102  mm  Embedded  Anchor 

the  progression  of  failure  for  the  127  mm  embedded  anchor.  Again  the  shaded  tension 
limit  initiated  at  the  concrete-adhesive  interface  and  propagated  toward  the  surface.  From 
the  boundaries  of  this  cone,  the  projected  failure  depth  was  62  mm  with  cone  radius  of 
1 1 7 mm.  Note  the  concrete  failure  cone  depth  is  also  decreasing  with  anchor  embedment 
depth.  The  gauss  point  at  which  the  element  stress  did  not  converge  was  62  mm  from  the 
top  of  the  mesh,  which  is  indicated  on  Figure  3.2 1 . The  point  is  inside  the  tension  zone. 

152  mm  embedment  anchors  Presented  in  Figure  3.22  is  the  simulation  of  the 
152.4  mm  embedded  anchor.  The  dots  characterize  either  tension  or  shear  limits.  The 
gauss  point  at  which  the  element  stress  did  not  converge  was  85  mm  from  the  top  of  the 
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Figure  3.19  Plan  View  of  Concrete  Cone  Failure  for  127  mm  Embedded  Anchor 


Figure  3.20  Cross-sectional  View  of  Concrete  Cone  Failure  for  127  mm  Embedded 

Anchor 
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Figure  3.21  Numerical  Simulation  of  Failure  for  127  mm  Embedded  Anchor 


mesh  and  the  dark  line  indicated  the  tension  failure  cone.  Also  displayed  in  the  hatched 
region  is  the  range  of  experimentally  observed  failure  surfaces. 

Table  3.3  shows  the  range  of  predicted  and  measured  concrete  failure  cones  for  all 
embedment  depths.  Also  shown  is  the  ratio  of  the  average  predicted  cone  depth  to  anchor 
embedment  depth  expressed  as  a percentage.  From  the  last  column  in  the  table  states 
that  the  percentage  of  the  overall  anchor  depth  associated  with  the  concrete  failure  cone 
decreases  with  increasing  anchor  length  with  the  exception  of  152  mm  embedded  anchor. 
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Figure  3.22  Numerical  Simulation  of  Failure  for  152  mm  Embedded  Anchor 
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Table  3.3  Measured  and  Predicted  Concrete  Failure  Cone  Dimensions 


Embedment 

Exp. 

Num. 

Exp. 

Num. 

havAf 

depth  (mm) 

r,/r2 

r 

h./hj 

h 

Exp. 

Num. 

(mm) 

(mm) 

(mm) 

(mm) 

(%) 

(%) 

76.2 

124/164 

164 

76/76 

73 

100 

96 

102 

73/177 

90 

40/72 

62 

55 

61 

127 

59/133 

117 

37/56 

62 

37 

49 

152 

87/140 

132 

41/76 

88 

38 

58 

Cook  (1991)  reported  a similar  increase  for  the  150  to  174  mm  embedded  anchors  versus 
the  125  mm  to  149  mm  embedded  anchors. 

3.5.5  Bond  Stress  Distribution  at  Adhesive/Concrete  Interface 

Shown  in  Figure  3.23  (a)  is  the  shear  stress  profiles  along  the  concrete-adhesive 
bond  interface  for  different  levels  of  axial  load  in  the  case  of  the  76  mm  embedded 
anchor.  Note  at  low  loads  (i.e.  16.7  kN)  the  shear  stress  is  high  at  the  top  and  then 
decreases  with  depth.  The  hyperbolic  tangent  shaped  shear  distribution  at  the  lower  loads 
is  the  linear  elastic  results,  since  no  elements  have  yielded.  This  corresponds  to  the 
elastic  solution  presented  by  Cook  et  al.,  (1993).  At  the  higher  loads  (33  kN  and  50  kN), 
both  the  epoxy  adhesive  and  the  concrete  begin  to  yield,  redistributing  the  load  (stress)  to 
the  bottom  of  the  adhesive  layer  (see  33.4  kN  load).  The  latter  occurs  whenever  the  yield 
surface  is  engaged,  i.e.  the  material  stiffness  decreases  over  it's  elastic  value  in  the 
yielded  area.  The  latter  explains  the  uniform  distribution  of  shear  along  the  epoxy- 
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concrete  interface  (8  MPa)  for  the  33.4  IcN  load.  Next,  as  the  load  approaches  failure,  the 
concrete  mass  near  the  surface  begins  to  dilate  (expand)  as  predicted  by  Figure  1.1  which 
increases  the  radial  stress  (also  mean  stress)  in  the  adhesive  layer  and  increases  its  failure 
shear  stress  (Figure  1.1)  as  shown  in  Figure  3.23  (a)  (see  33.4  kN  load). 

Presented  in  Figure  3.23  b,  c,  d,  e,  f,  g and  h is  the  developed  shear  stress  along 
the  102  mm,  127  mm,  152  mm,  178  mm,  254  mm,  356  mm  and  508  mm  embedded  16 
mm  diameter  anchors  with  different  levels  of  vertical  pullout  loads.  Again  quite  evident 
from  all  figures  is  the  increase  of  shear  stress  with  depth  until  a uniform  8 MPa  is  reached 
along  the  whole  anchor.  Then  the  materials  begin  to  dilate  (expand)  resulting  in  higher 
radial  stresses  which  increase  the  adhesive  layers  shear  capacity.  When  failure  initiates 
at  the  adhesive  / anchor  interface  and  propagates  to  form  the  concrete  cone,  the  radial 
normal  stress  then  rapidly  drops  in  the  adhesive  bond  region,  decreasing  its  shear 
capacity  and  the  whole  anchor  subsequently  fails.  Evident  from  Figure  3.23  b,  c,  d,  e,  f 
and  g the  failure  shear  stress  profiles  approach  a very  uniform  condition  (i.e.  constant 
shear  stress  vs  depth)  for  the  longer  anchor  embedment  lengths  below  the  failure  concrete 
cone  region.  From  Figure  3.23  h,  it  can  be  seen  that  the  anchor  has  failed  before  the  full 
mobilization  of  shear  stress  in  the  adhesive.  The  shear  stress  has  mobilized  only  up  to  390 
mm  depth  which  corresponds  to  h^^  / d of  24.5.  This  trend  was  also  shown  in  Figure  3.12 
where  the  shear  stress  was  almost  constant  up  to  h^^ /d  of  25.  From  these  findings  it  can  be 
concluded  that  the  optimum  h,f /d  is  about  25  beyond  which  the  bond  stress  will  not  fully 
redistribute  along  the  entire  length  of  the  anchor. 


Normalized  Depth 
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Shear  Stress  at 
Epoxy  / Concrete  Interface 


16.7  KN  — 33.4  KN  — 50.1  KN 


(a) 

Shear  Stress  at 
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Shear  Stress  at 


^ 16.7  KN  33.4  KN  50.1  KN  — 66.8  KN  83.5  KN  | 
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Shear  Stress  at 
Epoxy  / Concrete  Interface 


16.7  KN  33.4  KN  — 50.1  KN  — 66.8  KN  83.5  KN  — 100.2  KN 
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Normalized  Depth 
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Shear  Stress  at 
Epoxy  / Concrete  Interface 


16.7  KN  33.4  KN  ^ 50.1  KN  66.8  KN 
— 83.5KN  100.2  KN  — 116.9  KN 

S) 


Shear  Stress  at 
Epoxy  / Concrete  Interface 


16.7  KN  ^ 50.1  KN  --  116.9  KN-*-  164.8  KN 


(f) 
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Shear  Stress  at 
Epoxy  / Concrete  Interface 


-^33.4KN  83.46  KN— 180.8  KN^  215  KN 

(g) 


Shear  Stress  at 
Epoxy  / Concrete  Interface 


16.7  KN  — 123  KN  191 J 235.8  KN 


(h) 

Figure  3.23  Development  of  Shear  Stress  Along  Concrete-Adhesive  Layer  Interface  for 

Various  Embedment  Depths 
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3.5.6  Comparison  of  Predicted  and  Experimental  Failure  Loads 

Table  3.4  presents  the  measured  and  predicted  pullout  loads  for  the  76  mm,  102 
mm,  127  mm,  and  152  mm  adhesive-bonded  anchors.  The  column  with  the  heading  Exp. 
gives  the  experimentally  measured  values,  and  the  column  labeled  Num.  is  the  finite 
element  prediction  (maximum  error  1 1%).  The  column  identified  as  Eq.  3.7  is  the 
prediction  based  on  the  elastic  bond  stress  model  with  cone  depths,  hjo„j,  obtained  from 
experimental  data  (see  Table  3.3)  and  X obtained  from  Cook  (1991).  Note  the  major 
difference  (30%)  for  this  case  occurs  for  the  shallow  embedded  depths  were  the  concrete 
failure  cone  contribution  is  significant.  Recently,  Eligehausen  et  al  (1995)  has  modified 
the  full  embedment  anchor  pullout  force  given  in  Equation  3.1,  to 

{newtons)  (3.9) 

pullout  resistance  is  given  as 

=15.93  {newtons)  (3.10) 

The  column  identified  as  Eq.  3.10  is  the  prediction  based  on  Equation  3.10  with  the 
experimental  cone  depths  given  in  Table  3.3,  and  average  shear  stress,  Tq,  given  in  Table 
3.4.  From  the  comparison  of  measured  to  predicted,  it  is  evident  that  Equation  3.10 
results  in  too  high  (30%)  of  a pullout  resistance  at  the  lower  embedment  depths.  The 
final  column  is  identified  as  the  uniform  bond  model,  and  is  obtained  by  simply 
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multiplying  the  average  shear  (1 1 .5  MPa)  times  the  surface  area  of  the  anchor  hole  (7t 
hj^).  The  average  shear  stress  (1 1.5  MPa)  was  obtained  as  follows:  the  shear  stress  values 
at  each  gauss  point  along  the  interface  was  summed  up  and  divided  by  the  total  number  of 
gauss  points.  It  is  proposed  that  may  be  obtained  directly  from  a unconfined  pullout 
tests  on  anchors  with  different  embedment  depth  to  diameter  ratios  for  different  adhesive 
bonding  materials. 


Table  3.4  Measured  and  Predicted  Anchor  Pullout  forces  for  Different  Embedment 


Depths 


hcff 

V 

X 

Num. 

Exp. 

Eq.  3.7 

Eq.3.10 

Uniform 

(mm) 

(MPa) 

(MPa) 

(kN) 

(kN) 

(kN) 

(kN) 

(kN) 

76 

12.6 

14.5 

0.02 

52 

50 

33 

66 

52 

102 

11.65 

13.5 

0.02 

64 

72 

56 

73 

70 

127 

12.17 

14 

0.02 

82 

84 

77 

87 

87 

152 

11.71 

14 

0.02 

104 

105 

95 

107 

104 

178 

10.97 

13 

0.02 

111 

N.A 

N.A 

N.A 

121 

Note:  ^ Obtained  from  Figure  3.23 
N.A  Not  available 
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3.6  Summary 

The  following  conclusions  were  made  from  numerical  and  experimental  work: 

(1)  The  numerical  analysis  suggests  that  the  failure  initiates  at  the  concrete- 
adhesive  bond  interface  and  propagates  toward  the  siuface  to  form  the  shallow 
concrete  cone; 

(2)  The  limit  state  within  the  concrete  cone  is  tension  and  below  it  on  the 
concrete-adhesive  bond  interface  is  shear; 

(3)  Up  to  hjf /d  of  25,  the  ultimate  load  increases  proportionally  with  an  increase 
in  embedment  depth  and  ultimate  shear  stress  mobilized  at  the  interface  is  nearly 
constant.  Beyond  this  point  the  anchor  failed  before  full  mobilization  of  the  shear 
stress  in  the  adhesive; 

(4)  For  a given  embedment  depth  the  increase  in  ultimate  load  is  insignificant  for 
different  adhesive  thicknesses; 

(5)  For  a given  embedment  depth  the  increase  in  ultimate  load  is  insignificant  for 
different  concrete  strengths; 

(6)  For  a given  embedment  depth  the  ultimate  load  dropped  by  half  by  decreasing 
the  adhesive  shear  strength  by  half; 

7)  As  the  concrete  and  adhesive  layer  approach  failure,  they  dilate  (expand) 
resulting  in  an  increased  shear  capacity  of  the  adhesive  layer  which  is  lost  (no 
confinement)  once  the  failure  initiates  at  the  anchor  and  propogates  to  form  the 


concrete  cone;  consequently,  it  is  proposed  that  both  the  failure  cone  and 
concrete-adhesive  bond  fail  simultaneously; 

(8)  The  deeper  embedment  depth  (508  mm)  anchor  failed  before  full 
mobilization  of  the  shear  stress  in  the  adhesive.  The  uniform  shear  stress  was 
mobilized  only  up  to  390  mm  which  corresponds  to  h^f /d  of  25.  This  clearly 
shows  that  uniform  bond  stress  can  be  applied  up  to  h^f /d  of  25; 

(9)  In  practice  h^f /d  is  limited  to  the  order  of  20-25,  at  this  level  full  redistribution 
of  shear  stress  is  taking  place  as  shown  in  Figures  3.23  and  the  uniform  bond 
stress  model  can  be  used  to  predict  the  capacity  of  the  adhesive-bonded  anchor. 

In  order  to  test  the  uniform  bond  model  for  larger  h^f /d  a bond  slip  element  may 
be  needed  to  predict  the  capacity  accurately,  this  is  discussed  in  the  next  chapter. 


CHAPTER  4 

DEVELOPMENT  OF  AXISYMMETRIC  QUADRATIC  INTERFACE  ELEMENT  OF 

PENALTY  TYPE 


4.1  Introduction 

In  recent  years,  there  has  been  a considerable  interest  in  application  of  the  finite 
element  method  to  discontinuous  systems,  especially  in  the  field  of  rock  mechanics  where 
the  discontinuities  consist  of  rock  joints,  faults  and  various  geological  interfaces.  Other 
engineering  problems  where  interface  discontinuities  play  an  important  role  are  joint 
connections  in  composite  structures,  analysis  of  fuel  element  assemblies  in  nuclear 
reactors,  and  machine-part  interfaces,  such  as  gears  and  rollers. 

Another  important  class  of  discontinuity  problems  that  civil  engineers  are  faced 
with  are  the  inelastic,  nonlinear  behavior  found  in  soil-structure  interaction  (such  as  piers, 
footings,  culverts),  masonry  structures,  load-carrying  mechanism  between  concrete  and 
reinforcement  in  the  longitudinal  direction  of  the  reinforcing  bars,  the  bond  between 
steel/adhesive  and  adhesive/concrete  in  chemically  bonded  anchors,  etc.  Such  interfaces 
constitute  “planes  of  weakness”  and  the  ability  to  describe  their  behavior  with  reasonable 
accuracy  is  essential. 

By  using  regular  finite  element  methods  to  analyze  the  above  mentioned  problems 
without  any  special  provisions  to  model  such  discontinuities,  then  the  elements  would  be 
connected  to  each  other  at  their  respective  nodal  points.  Consequently,  two  adjacent 
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elements  would  have  identical  displacements  at  the  shared  nodes  and  a smooth  stress 
field  would  develop.  However,  two  dissimilar  materials  would  generally  have  a jump 
discontinuity  in  the  stress  field  across  their  boundaries.  The  latter  could  only  occur,  if  the 
two  materials  had  different  strains,  which  would  require  relative  nodal  displacements. 
Consequently,  if  bond-slip  is  to  be  taken  into  account,  special  interface  elements  with 
zero  thickness  must  be  used  for  robust  idealization  of  the  interfaces  between  the  two 
dissimilar  materials. 

The  physical  behavior  of  an  interface  separating  two  continua  involves  tangential 
slip  initiation  when  the  tangential  component  of  the  interface-stress  vector  exceeds  the 
“shear  strength”  of  the  discontinuity.  Subsequently  relative  slip  according  to  some 
“sliding  friction”  law,  or  separation,  or  both  will  occur  if  the  normal  component  of  the 
interface-stress  vector  exceeds  the  “tensile  strength”  of  the  discontinuity.  Slip, 
separation,  and  recontact  may  all  occur  during  a given  stress-time  history. 

The  interface  problem  is  complex  and  inherently  nonlinear.  From  a mathematical 
standpoint,  analytical  closed  form  solutions  are  possible  only  for  a limited  class  of 
idealized  interface  problems  (Duvant  and  Lions  1976).  As  such,  there  has  been  a growing 
interest  in  numerical  solution  capabilities  for  such  problems. 

Various  numerical  methods  have  been  proposed  in  the  literature  to  solve  interface 
problems.  Researchers  like  Zienkiewicz  et  al.  (1970)  and  Pande  and  Sharma  (1979)  have 
treated  interface  as  a quasi-continuum  of  small  thickness  i.e.  by  using  continuum  finite 
elements  which  contain  planes  of  weakness.  Others  like  Francavilla  and  Zienkiewicz 
(1975),  Sachdeva  and  Ramakrishnan  (1981),  Chan  and  Tuba  (1971)  and  Hughes  et  al. 
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(1976)  have  modeled  discontinuous  behavior  using  constraint  equations  or  by  connecting 
the  elements  with  each  other  by  discrete  springs  also  known  as  stiffiiess  method 
(Hermann,  1978;  Frank  et  al.,  1982;  Goodman  et  al.,  1968;  Ghaboussi  et  al.,  1973;  and 
Schafer,  1995).  The  latter  work  is  of  particular  interest,  and  these  are  known  as  zero 
thickness  elements. 

The  method  of  stiffiiess  is  basically  the  simple  concept  of  using  “bar”  elements 
(or  directionally  stiff  elements)  across  the  interface  in  both  the  normal  and  tangent 
direction.  For  example,  to  model  frictionless  slippage  across  an  interface,  the  normal 
stiffness  would  be  specified  arbitrarily  large  to  force  near  compatibility  of  normal 
displacements,  while  the  tangent  stiffness  would  be  specified  extremely  small  (or  zero)  to 
allow  independent  movement  in  the  tangent  direction.  This  large  stiffness  is  required  to 
prevent  nodes  on  either  side  of  the  zero-width  interface  from  penetrating  each  other  under 
compressive  loading.  On  the  other  hand  if  a normal  stiffness  is  selected  too  large  with 
respect  to  the  computer  word  length,  significant  digits  of  the  relative  displacement 
become  truncated,  and  the  calculation  of  interface  forces  may  be  in  error.  If  the  normal 
stiffiiess  is  selected  too  low,  significant  penetration  may  occur,  and  the  kinematics  may  be 
in  error  (body  penetration).  Also  Gens  et  al.  (1989)  and  Schellekens  and  De  Borst  (1993) 
have  cited  other  inconsistencies:  spurious  stress  oscillations  attributed  due  to 
inappropriate  quadrature  schemes,  lack  of  accuracy  due  to  excessively  large  stiffness 
parameters  and  inaccurate  interface  stress  predictions  due  to  insufficient  mesh  fineness. 
Kaliakin  et  al.  (1995)  in  a recent  article  have  studied  the  linear  models  of  Goodman  et  al. 
(1968)  and  Herman  (1978)  in  detail.  The  authors  have  found  that  the  inconsistencies 
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were  not  because  of  the  factors  mentioned  above,  but  were  due  to  certain  fundamental 
kinematic  deficiencies  associated  with  the  zero-thickness  interface  elements  themselves. 
On  eliminating  the  fundamental  deficiencies  , the  zero-thickness  elements  were  taken 
through  several  examples.  The  improved  zero-thickness  element  showed  superior 
behavior  compared  to  the  Goodman  and  Herman  elements. 

The  method  of  constraints  is  a more  accurate  than  the  stiffness  method.  Here  an 
interface  discontinuity  in  a discretized  model  of  a structure  is  represented  by  a sequence 
of  double  nodes,  one  on  each  side  of  the  interface.  The  interconnection  between  double 
nodes  is  introduced  through  kinematics  and  the  desired  solution  is  obtained  by  modifying 
the  global  stiffness  equations  in  such  a manner  that  all  the  interface  conditions, 
compatibility  and  friction  law,  are  satisfied.  These  additional  constraints  makes  it  more 
expensive  computational  wise,  especially  when  solving  large  problem. 

In  this  chapter,  an  improved  quadratic  axisymmetric  zero-thickness  element  will 
be  derived  using  the  method  of  stiffhess  and  implement  into  the  general  purpose  non- 
linear finite  element  program  PLASFEM.  This  element  will  be  free  from  some  of  the 
inconsistencies  mentioned  earlier. 

4.2  Derivation  of  Element  Stiffhess  Matrix 

When  using  quadratic  continuum  elements  in  two-dimensional  analyses,  nodal 
compatibility  requires  the  use  of  quadratic  zero-thickness  interface  elements.  This 
element  possess  six  nodes  and  twelve  displacement  degrees  of  freedom  Figure  4.1.  At 
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each  nodal  location  along  the  interface  a pair  of  nodes  is  placed  at  the  same  initial 
geometric  location;  the  thickness  of  the  element  is  thus  initially  zero. 

The  vector  of  absolute  nodal  displacement  degrees  of  freedom  with  respect  to  the 
local  coordinate  system  (s-n ) is  given  by 

{u}  = {ui  Vi  U2  V2  U3  V3  U4  V4  U5  V5  U6  V6  (4.1) 

where  the  superscript  T denotes  the  operation  of  matrix  transposition. 


Z 


Figure  4.1  Six-Node  Axisymmetric  Zero-Thickness  Element  of  Goodman  et  al.  (1968) 
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Associated  with  the  above  displacement  components  is  a corresponding  vector  of 
nodal  forces 

{f'}  = {qi  q2  q3  q4  qs  qe  q?  q»  q9  qio  qii  qi2  (4.2) 

The  vector  of  relative  displacements  is  defined  as 

= (4.3) 


where  Ws  and  w„  represent  tangential  and  normal  displacements  respectively  along  the 
interface.  The  matrix  relating  relative  and  absolute  displacement  is  given  by 


M= 


-N, 


0 -N, 


0 -N, 


0 -N,  0 -N, 


0 Nj  0 0 0 

- A^3  0 A^3  0 0 


where  Ni , N2  and  N3  represent  standard  quadratic  interpolation  functions  associated  with 


the  ^ - coordinate  system  and  they  are  given  as 


(1-#  ) 

iV2=(l-^^  (4.4) 

(l  + « ) 


N4 , N5  and  Ne  are  same  as  N3 , N2  and  Ni. 

The  relative  displacements  are  related  to  forces  acting  along  the  interface  through  the 
following  general  constitutive  relation 


'dn 

1 

A2 

^22. 

Kj 

(4.5) 
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M=[D]{w} 


(4.6) 


where  CTs  and  On  represent  the  tangential  and  normal  forces  per  unit  length,  respectively, 

acting  along  the  interface.  The  constitutive  parameters  dn  and  d:2  represent  the 
tangential  and  normal  stiffness  per  unit  length,  respectively,  along  the  interface.  di  i and 
d22  can  also  be  thought  of  as  penalty  numbers  along  tangential  and  normal  directions. 
The  presence  of  the  coupling  parameters  di2  allows  for  the  simulation  of  dilatational 
phenomena. 


If  the  mid-side  nodes  are  equidistant  from  the  comer  nodes  and  if  the  length  of  the 
element  is  denoted  by  L,  the  radius  r and  the  | reduces  to 


The  elemental  stiffness  matrix,  is  computed  from 


(4.7) 


(4.8) 


(4.9) 


After  performing  the  multiplication  of  [N*]^  [D]  [N*] , the  final  product  looks  like 
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(4.10) 
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iV13 
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N22  (J22 
0 

N33  d^2 


-A13 
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-N12 

0 

-Ml  d,, 

0 
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-N\3  d22 

0 

-N\2  d22 

0 

- m 1 d22 

-N23 

0 

-N22  c/„ 

0 

-M2 

0 

0 

- N23  d22 

0 

-N22  d22 

0 

-M2  d22 

-n33 

0 

-N32  d^^ 

0 

-Ml  d„ 

0 

0 

- N33  d22 

0 

-N32  d22 

0 

-Ml  d22 

0 

M2  0 

Ml  c/„ 

0 

A^33  d22 

0 N32  d22 

0 

A31  d22 

N22  0 

Ml 

0 

N22  d22 

0 

A^21  d22 

Sym. 

Ml 

0 

N\  1 d22 

By  looking  at  the  various  terms  in  the  above  matrices,  (A , B,  and  C)  there  are  six 
different  integrals  which  have  to  be  evaluated.  These  are  N1 1,  (which  means  N1 
multiplied  by  Nl),  N12,  N13,  N22,  N23,  N33.  An  explicit  integration  is  performed.  A 


A B 

y C 

The  matrices  A,  B,  and  C are  defined  as 


0 

N\2 

0 

M3  d, 

N\  1 d22 

0 

N\2  d22 

0 

N22  c/„ 

0 

N23  d, 

N22  d22 

0 

Sym. 

N33  d, 
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sample  calculation  for  one  of  them  is  given  and  for  the  rest  of  the  integrals,  only  the  final 
results  are  given.  The  integrated  values  are  then  multiplied  with  the  associated  normal 
and  tangential  stiffness  to  get  the  joint  element  stiffness  per  unit  length  which  is  in  the 
element  local  coordinate  system. 

Ml  = A^l*^l  = ;ri:j  [-|(l-^)f (F + | 

N11  = ^(8F-3L) 

N22  = ^(32r) 

N33  = ^(8r+3L)  (4.11) 

N12  = ^(4r-2L) 
n L 

iV13  = — (-2r) 

N23  = ^{4r^2L) 

The  last  step  in  the  joint  stiffness  derivation  is  to  place  the  element  in  a coordinate 
system  for  the  entire  structure.  Adopting  a global  coordinate  system  (r,  z),  as  drawn  in 
Figure  4.1,  the  two  systems  are  related  through 


cos^ 

-sin^ 


sin^ 

cos^ 


(4.12) 
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It  follows  that  the  element  stiffiiess  matrix  and  force  vector  with  respect  to  global 
coordinates  are  computed  in  the  usual  manner,  that  is 


[/‘1=[«T[/1 


where 


[r]= 


'[R]  0 0 0 0 o' 

0 [/?]  0 0 0 0 

0 0 [/?]  0 0 0 

0 0 0 [/?]  0 0 

0 0 0 0 [/?]  0 

0 0 0 0 0 [/?] 


(4.13) 

(4.14) 


(4.15) 


4.3  Behavior  of  a Single  Quadratic  Axisvmmetric  Zero-Thickness  Element 

In  order  to  gain  insight  into  the  behavior  of  the  zero  thickness  element,  consider 
the  single  element  example  shown  in  Figure  4.2.  In  this  case,  the  local  element  axis  is 
aligned  with  the  global  axis  for  simplicity. 

Accounting  for  the  homogenous  displacement  specifications  at  nodes  1, 2 and  3 


the  element  equation  reduces  to 


102 


■ IdW 

0 

AdW 

0 

-d\\ 

0 

«4 

^7 

0 

Idll 

0 

Adll 

0 

-dll 

V4 

7d} 

Ad\\ 

0 

\6d\\ 

0 

0 

0 

^9 

0 

Adll 

0 

\6dll 

0 

0 

V5 

* = ’ 

9l0 

-~d\\ 

0 

0 

0 

d\\ 

0 

^11 

0 

-dll 

0 

0 

0 

dll  _ 

.V'6. 

.^12. 

(4.16) 


For  a special  case  of  q;  = P , q9  = qi  i = 0 and  qg  = qi2  = - Q , qio  = - 4Q,  the  nodal 
displacements  are  given  as 


L^6 


f ^ 

d\\ 

-6Q 


TdJ 


(4.17) 


dll 
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Figure  4.2  Example  Involving  Single  Original  Zero-Thickness  Element 
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The  quantities  U4,  U5  and  U6  represent  the  “pre-slip”  displacements;  i.e  those 
occuring  prior  to  the  mobilization  of  the  frictional  resistance  along  the  interface.  Even  if, 
through  the  choice  of  large  values  for  di  1,  these  displacements  are  minimized,  the  fact 
remains  that  a kinematic  inconsistency  is  associated  with  the  element.  That  is,  node  5 
displaces  in  a direction  opposite  to  that  associated  with  the  applied  load  P at  node  4 
(Figure  4.2).  The  normal  response  is  seen  to  be  free  of  any  kinematic  difficulties. 

4.4  Development  of  Improved  Quadratic  Axisvmmetric  Interface  Element 

In  the  light  of  the  above  deficiencies,  it  is  evident  that  an  improved  zero-thickness 
interface  element  is  desirable.  The  development  of  such  an  element  is  now  presented. 

Consider  a “macro-element”  consisting  of  two  6 node  elements  as  shown  in 
Figure  4.3.  The  thickness  is  initially  zero.  The  pairs  of  nodes  1 and  6, 2 and  5,  3 and  4,  7 
and  8,  and  9 and  10  are  coincident  before  deformation.  There  are  no  external  nodal 
forces  acting  at  nodes  7 through  10.  As  given  in  Section  4.2  for  each  element  the  stress- 
displacement  is  the  same.  For  each  element  the  resulting  stiffness  matrix  with  respect  to 
the  local  element  coordinates  is  thus  given  by  Equation  4.13  with  proper  L’s  and  r’s  for 
each  element. 

The  nodal  displacement  and  nodal  force  vectors  for  each  element  are: 

For  element  1 

{u}i  = {ui  Vi  U7  V7  U2  V2  U5  V5  Ug  Vg  U6  V6  (4.18) 


92  qi3  qi4  qs  q4  q9  qio  qis  qi6  qii  qi2  (4.19) 
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For  element  2 


{u}2  = { U2  V2  U9  V9  U3  V3  U4  V4  UiQ  Viq  U5  V5 


Q2  qi3  qi4  q3  q4  q9  qio  qis  qi6  qii  qi2 


where  the  superscript  T denotes  the  operation  of  matrix  transposition. 


r 


(4.20) 

(4.21) 


Figure  4.3  “Macro-Element”  Consisting  of  Two  Axisymmetric  Zero-Thickness  Elements 
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To  form  the  element  stiffness  equations  of  the  two  elements  in  local  coordinates, 
the  stiffiiesses  of  the  two  elements  are  added  using  the  standard  direct  stiffiiess  approach. 
The  stiffness  of  the  system  thus  have  the  following  form 

\KY{uV={fV  (4.22) 

where  superscript  A means  assembled  matrices.  The  right  had  side  of  the  system 
equations  has  also  been  assembled  from  the  right  had  side  vector  {f  *}’s  of  each  element. 
The  forces  at  the  internal  nodes  are  zero.  The  assembled  force  is  given  below, 

{f  = { qi,  q2,  qs,  q4,  qs,  q6,  q?,  q8,  q9,  qio,  qii,  qi2,  0,  0,  0,  0,  0,  0,  0,  0}  (4.23) 

If  the  mid  side  nodal  displacements  for  each  element  are  expressed  in  terms  of  the 
comer  nodal  displacements  of  the  respective  element  we  get  the  following  relationships: 
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V9-V,o  = 


0 0 0 -- 


,A^12 


0 - 


2N23 

2^2' 


0 


2NI2 

2^2^ 


0 0 


{«} 


where  {u}  is  the  vector  of  1 through  12  nodal  displacements,  which  are  same  as 
defined  in  Equation  4.1 


{u}  = { Ui  Vi  U2  V2  U3  V3  U4  V4  U5  V5  U6  V6 


(4.24) 


Next  the  element  equations  shall  be  condensed,  reducing  the  nodal  displacement 
vector  {u}''  into  {u}.  Expanding  the  twelve  equations  corresponding  to  {u}  and 
rewriting  the  terms  with  {u}  into  a matrix  equation  give 

[/:](«) +{G)  = {/1  (4.25) 

where  {G}  = { gu,  gin,g2s,  g2n,  g3s,  g3n,  g4s,  g4n,  gss,  gsn,  g6s,  g6n}  is  a vector  Containing 
terms  associated  with  U7,  V7,  ug,  vs,  U9,  V9,  uio  and  viq.  We  can  easily  find  the  expression 
for  {G}  by  expanding  the  corresponding  equations  in  the  same  rows  in  Equation  4.22  and 
Equation  4.25  and  comparing  their  terms: 
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Sis  ~ 
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The  terms  in  G matrix  associated  with  the  normal  direction  are 
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Substituting  these  terms  into  the  G matrix,  and  summing  with  the  [K]  matrix  the  final 
condensed  stiffness  matrix  in  the  local  coordinate  is 
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(4.26) 


The  transformation  of  element  equations  to  global  system  is  realized  in  the  manner 
described  by  Equation  4.15. 

4.4.1  Interface  Stresses 


The  interface  tangential  and  normal  forces  per  unit  length  are  defined  on  the 
following  manner 


<j,  = 


^7  +^9+^11 

rL 


(4.27) 


cr.  = 


^8  ^10  ^12 
rL 


(4.28) 


where  equation  4.26  has  been  used  to  substitute  for  q7  through  qi2. 


4.4.2  Predictive  Capabilities  of  Improved  Quadratic  Axisvmmetric  Zero-Thickness 
Element 

The  single  element  example  shown  in  Figure  4.2  is  again  considered.  A very 
large  radius  of  1000  units  was  chosen  to  simulate  plane  strain  condition.  Incorporating 
the  homogeneous  nodal  specifications  and  solving  the  set  of  equations  from  4.26  gives 
horizontal  and  vertical  displacements  as 
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Contrary  to  the  tangential  response  of  the  original  quadratic  axisymmetric  zero-thickness 
element  (Equation  4.17)  where  the  mid-side  node  displaced  in  a direction  opposite  to 
applied  load,  the  tangential  displacements  associated  with  the  improved  quadratic 
axisymmetric  zero-thickness  are  seen  to  be  consistent.  As  such,  the  kinematic  deficiency 
noted  for  the  original  element  has  been  removed  from  the  formulation. 


4.5  Example  Involving  Simulated  Pull-Out  Experiment 

In  order  to  further  assess  the  behavior  of  the  improved  quadratic  axisymmetric 
element  and  to  compare  with  the  Kaliakin  improved  two  node  plane  strain  element,  an 
example  involving  the  simulated  pull-out  of  an  inclusion  from  a surrounding  continuum 
is  next  analyzed.  A typical  finite  element  mesh  used  and  the  boundary  conditions  are 
shown  in  Figure  4.4.  The  continuum  is  modeled  by  eight  noded  axisymmetric  elements 
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Figure  4.4  Simulated  Pull-out  Example  Involving  Zero-Thickness  Interface  Elements 

numbered  1 through  48.  Elements  49  to  72  are  interface  elements,  and  the  inclusion  is 
represented  by  twelve  eight  noded  axisymmetric  elements  (with  an  aspect  ratio  of  2), 
Both  the  continuum  and  the  inclusion  are  idealized  as  isotropic,  linear  elastic  materials. 
The  elastic  modulus  and  poisson’s  ratio  for  the  continuum  are  2.0  x 10^  and  0.25 
respectively.  The  self  weight  of  the  continuum  material  is  considered  to  be  negligible, 
and  a state  of  plane  is  assumed,  for  this  purpose  the  radius  is  assumed  of  the  order  of 
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10000  units.  For  the  inclusion,  the  elastic  modulus  and  thickness  are  taken  equal  to  2.0  x 
lO'*  and  0.25  respectively.  The  shear  strength  of  the  interface  is  assumed  to  be  governed 
by  Mohr-Coulomb  criterion  with  a friction  coefficient  equal  to  0.50  and  cohesion 
intercept  equal  to  zero.  A normal  distributed  compressive  load  equal  to  10.0  is  applied  in 
the  first  solution  step  and  is  maintained  constant  throughout  the  remainder  of  the  loading 
history.  The  horizontal  load  is  applied  monotonically  from  a value  of  zero  at  the  first 
solution  step  to  a value  of  50.0  at  the  25  step. 

The  normal  response  of  the  interface  elements  is  first  assessed.  Since  in  the  first 
solution  step  only  the  uniform  normal  distributed  load  is  applied  to  the  mesh,  the  normal 
respx)nse  can  easily  be  assessed  from  just  this  one  step.  The  analysis  was  carried  out  for 
different  values  of  d2i  such  as  10*  and  lO'^.  Even  for  these  high  values  622  the  normal 
stresses  are  exactly  transmitted  across  the  interface  and  the  relative  normal  displacements 
across  the  interface  are  effectively  zero  of  the  order  of  1 .0  x 10’*. 

The  combined  normal  and  tangential  response  of  the  element  is  considered  next. 
The  values  of  dn  = 1.0  x 10^  and  ^22  - 4.0  x 10'°  are  assumed.  Kaliakin  reported  results 
for  6,  12, 24  and  48  elements  per  side.  It  is  found  that  results  were  identical  for  meshes 
containing  12, 24  and  48  elements  per  side.  Here  the  mesh  selected  contains  12  interface 
elements  per  side.  The  results  are  presented  in  Figure  4.5.  From  the  figure  it  is  evident 
that  the  results  are  practically  identical  to  those  reported  by  Kaliakin.  The  solution 
converged  in  2-3  iterations. 
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In  summary,  this  example  has  shown  that  the  normal  response  is  quite  similar  to 
that  associated  with  the  Goodman  element  and  is  free  from  the  force  oscillations 
associated  with  the  latter  element. 


Distance  Along  Inclusion 


Figure  4.5  Distribution  of  Total  Tangential  Force  Along  Inclusion  at  Load  of  30.0 


CHAPTER  5 
CONCLUSIONS 

5.1  Summary 

Development  of  a robust,  economical  and  unconditionally  stable  nonlinear  finite 
element  program  PLASFEM  is  an  ongoing  project  in  Geotechnical  Engineering, 
Department  of  Civil  Engineering,  University  of  Florida.  The  program  is  PC  based  and 
can  handle  from  2-D  and  3-D  problems.  The  program  also  includes  various  constitutive 
models.  These  models  can  be  used  to  model  a number  of  material  nonlinear  responses. 
The  program  was  deficient  in  modeling  time  dependent  behaviors  and  interface  problems. 
In  order  to  handle  the  former,  a viscoplastic  model  of  Perzyna  type  was  formulated  and 
implemented  for  the  Sandler  and  Dimaggio  cap  model.  The  program  is  unconditionally 
stable  for  any  load  step  size.  The  above  mentioned  features  were  achieved  by  directly 
evaluating  the  consistent  tangent  operator  at  each  gauss  point  per  iteration  without  any 
costly  material  stiffness  inversions.  This  ensures  a quadratic  rate  of  convergence  for  the 
Newton  solution  procedure  (Equation  1 .6).  Since  the  yield  surface  is  composed  by  a 
tension  cutoff,  failure,  a hardening  elliptic  cap  and  comers  (tension  and  compression), 
four  different  operators  were  developed  none  of  these  required  any  material  or  stif&iess 
inversions  and  were  easily  implemented  into  the  finite  element  code  PLASFEM.  The 
converged  stresses  due  to  any  load  increment  were  found  using  an  implicit  backward 
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Euler  formulation,  usually  referred  as  the  closest  point  projection.  The  updated  stresses 
determined  were  in  terms  of  stress  invariants;  I,  (P‘  invariant  of  stress)  and  V (2™* 
invariant  of  deviatoric  stress). 

In  order  to  model  discontinuity  problems  faced  by  civil  engineers,  such  as  soil- 
structure  interaction  (such  as  piers,  footings,  culverts),  masonry  structures,  bond  slip 
between  concrete  and  reinforcement  in  the  longitudinal  direction  of  the  reinforcing  bars, 
the  bond  between  steel/adhesive  & adhesive/concrete  in  chemically  bonded  anchors,  etc., 
an  interface  element  is  required  to  model  two  dissimilar  materials.  An  axisymmetric 
quadratic  interface  element  of  penalty  type  was  formulated  and  implemented  into 
PLASFEM  to  handle  the  above  mentioned  problems. 

An  extensive  study  of  behavior  of  single  adhesive-bonded  anchors  was  made 
using  the  viscoplastic  cap  model. 


5.2  Conclusions 

The  viscoplastic  cap  model  was  used  to  model  several  examples.  The  first  in  the 
series  was  McCormick  Ranch  Sand  example.  This  problem  was  used  to  study  the  effect 
of  different  fluidity  parameters  on  the  stress  response.  For  large  values  of  fluidity 
parameters  the  response  is  close  to  inviscid  (elasto-plastic)  throughout  the  history  and  for 
small  values  (close  to  zero),  the  response  is  close  to  pure  elastic.  The  next  example 
which  is  worth  noting  is  the  variable  loading  rates  on  sand  in  uniaxial  strain.  Upon 
examining  Figure  3.9  it  is  evident  that  the  viscoplastic  cap  model  accurately  reflects  the 
test  data.  For  the  slow  loading  test  the  cap  surface  continually  moves  with  the  stress  state 
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producing  a stiffening  stress/strain  response  as  shown.  On  the  other  hand,  for  the  fast 
loading  rate  the  cap  surface  lags  behind  the  loading  stress  state  producing  an  apparent 
softening  stress/strain  response  which  is  truly  a time-dependent  effect.  The  correlation 
between  model's  performance  and  the  observed  performance  is  truly  remarkable. 

Next  the  model  was  used  to  study  the  inviscid  behavior  of  single  adhesive  anchors 
under  tensile  load.  The  experimental  effort  involved  18  unconfined  pullout  tests  having 
four  different  embedment  depths  (76  mm,  102  mm,  127  mm,  and  152  mm)  on  one 
adhesive  bonding  agent  (an  epoxy  amine)  and  one  diameter  threaded  rod  (16  mm).  The 
following  important  conclusions  were  arrived  at  based  on  both  the  numerical  and 
experimental  work: 

1)  The  numerical  analysis  suggests  that  the  concrete  failure  cone  initiates  at  the 
concrete-adhesive  bond  interface  and  propagates  toward  the  surface; 

2)  The  observed  failure  cones  for  each  of  the  tested  embedments  compared  quite 
favorably  with  the  numerical  prediction; 

3)  The  limit  state  within  the  concrete  cone  is  tension  and  below  it  on  the  concrete- 
adhesive  bond  interface  is  shear; 

4)  As  the  concrete  and  adhesive  layer  approach  failure,  they  dilate  (expand) 
resulting  in  an  increased  shear  capacity  of  the  adhesive  layer  which  is  lost  (no 
confinement)  once  the  failure  initiates  at  the  anchor  and  propagates  to  form  the 
concrete  cone;  consequently,  it  is  proposed  that  both  the  failure  cone  and 
concrete-adhesive  bond  fail  simultaneously; 
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5)  Due  to  (4),  the  concrete-adhesive  bond  interface  strength  in  adhesive-bonded 
anchors  is  not  adequately  predicted  by  a direct  shear  test,  (i.e.  an  epoxied  plate 
pulled  along  a concrete  base).  Bond  strength  should  be  determined  from  testing  of 
chemically  bonded  anchors; 

6)  As  the  embedment  depth  increases,  the  shear  stress  on  the  concrete-adhesive 
bond  interface  becomes  more  uniform; 

7)  The  concrete  strengths  did  not  significantly  influence  the  ultimate  load  carrying 
capacity  of  the  anchor; 

8)  By  reducing  the  adhesive  shear  stress  the  capacity  of  the  anchor  dropped 
proportionately; 

9)  Increases  in  adhesive  thickness  increased  the  ultimate  load  by  only  about  5-10 
% for  a particular  diameter  and  embedment  depth; 

10)  Shear  stress  mobilized  at  the  interface  is  nearly  constant  at  about  13  MPa  up 
to  hjf  /d  of  25  and  then  drops; 

11)  Ultimate  load  increases  proportionally  with  an  increase  in  embedment  depth 
up  to  hjf/d  of  25; 

12)  From  (6),  and  experimental  data,  it  is  proposed  that  the  use  of  a uniform 
stress  model,  x,^g7tdo  heir » where  h^ff  is  embedment  depth  of  the  anchor,  and  T„g  is 
obtained  from  experimental  tests,  will  result  in  less  than  4%  error  for  the  4 
embedment  depths  investigated. 

Lastly  investigation  of  quadratic  axisymmetric  zero-thickness  interface  element 
similar  to  Goodman  et.  al.  (1968),  revealed  a fundamental  kinematic  deficiency  inherent 
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in  its  response.  In  past  investigations,  such  oscillations  were  attributed  to  the  particular 
numerical  integration  scheme  used  in  computing  the  element  stiffiiess.  However  it  is 
shown  here,  the  aforementioned  force  oscillations  can  be  removed  not  by  altering  the 
numerical  integration  scheme  but  by  removing  the  kinematic  deficiency  that  plagues  the 
element  through  condensation.  Based  on  the  above  observations,  an  improved  quadratic 
axisymmetric  zero-thickness  interface  element  was  developed.  This  element  possess  the 
normal  response  of  the  Goodman  model  but  eliminates  above  mentioned  kinematic 
deficiencies. 

The  improved  quadratic  axisymmetric  zero-thickness  element  was  verified 
through  pull-out  simulation  of  an  inclusion  from  a surrounding  continuum.  The  normal 
response  of  improved  element  was  shown  to  be  quite  similar  to  that  associated  with  the 
Goodman  model,  but  is  free  from  the  associated  pathologies.  The  improved  element  was 
also  shown  to  be  fi-ee  of  the  force  oscillations  that  Goodman  element  possessed.  Finally 
the  convergence  on  the  tangential  force  distributions  was  found  to  be  superior  to  the 
Goodman  element. 

In  short,  a robust  nonlinear  finite  element  program  which  is  under  development 
can  model  time  dependent  (creep),  time  independent  problems  and  discontinuity 
problems  faced  by  civil  engineers. 

5.3  Recommendation  for  Future  Work 

Since  the  program  PLASFEM  is  a general  purpose  program  it  can  be  used  to  study 
the  behavior  of  group  anchors,  edge  effects  and  creep  of  adhesive-bonded  anchors.  There 
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is  potential  for  adding  kinematic  and  isotropic  hardening  into  PLASFEM.  These  models 
can  be  used  to  model  cyclic  loading  such  as  earthquake  loading. 
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